Hello,
i need to integrate a very complicated function. I tried with scipy but the numerical integration is incorrect. So, i decided to use ROOT for numerical integration.
The function is a two-dimensional function, very long function! I decided to use python because i need to operate with matrices and trigonometric functions, if you suppose that it is better to use C++ please tell me.
The main problem is that i can’t implement my function inside a TF2. A very simple example:
import ROOT
import numpy as np
def myFunc(x,y):
return y*x
f = ROOT.TF2("f","myFunc(x,y)",0,10**8,0,2*ROOT.TMath.Pi())
wf1 = ROOT.Math.WrappedMultiTF1(f)
ig = ROOT.Math.AdaptiveIntegratorMultiDim()
ig.SetFunction(wf1)
ig.SetRelTolerance(0.001)
xmin = ROOT.std.vector('double')(2)
xmax = ROOT.std.vector('double')(2)
xmin[0] = 0
xmin[1] = 0
xmax[0] = 10**8
xmax[1] = ROOT.TMath.PiOver2()
print ig.Integral(np.array(xmin),np.array(xmax))
print f.Integral(0,10**8,0,ROOT.TMath.PiOver2())
This gives the error:
Error in <TFormula::Eval>: Formula is invalid and not ready to execute. myFunc is unknown.
The complete code is:
import ROOT
import numpy as np
from math import sin, cos
from sympy import simplify
from scipy import integrate
pi = ROOT.TMath.Pi()
e = ROOT.TMath.E()
mu0 = 4*pi*10**(-7)
e0 = 8.854*10**(-12)
c = 299792458
h = 1.05457*10**(-34)
wr = 1
wk, wm, wkr, wpl, gk, gm, gr, al, A = 0.7*wr, wr, wr, 0, 0.05*wr, 0.005*wr, 0.005*wr, 0.001, 0.001
gpl, we = 0.005*wpl, wkr
kk = lambda x,t: 1j*wk*(x*cos(t))/(-(x*cos(t))**2-wkr**2-gk*(x*cos(t)))
mu = lambda x,t: 1 + al + A*(x*cos(t))**2/(-(x*cos(t))**2-wm**2-gm*(x*cos(t)))
eps = lambda x,t: 1 - wpl**2/(-(x*cos(t))**2-gpl*(x*cos(t))) - we**2/(-(x*cos(t))**2-wr**2-gr*(x*cos(t)))
eta0 = np.sqrt(mu0/e0)
eta = lambda x,t: eta0*np.sqrt(mu(x,t)/eps(x,t))
n1 = lambda x,t: np.sqrt(eps(x,t)*mu(x,t)) + kk(x,t)
n2 = lambda x,t: np.sqrt(eps(x,t)*mu(x,t)) - kk(x,t)
gamm1 = lambda x,t: simplify((eta0**2 + eta(x,t)**2)/(2*eta(x,t)*eta0))
gamm2 = lambda x,t: simplify((eta0**2 - eta(x,t)**2)/(2*eta(x,t)*eta0))
chi1 = lambda x,t: simplify(np.sqrt((x*sin(t))**2+n1(x,t)*(x*cos(t))**2)/(n1(x,t)*x))
chi2 = lambda x,t: simplify(np.sqrt((x*sin(t))**2+n2(x,t)*(x*cos(t))**2)/(n2(x,t)*x))
delta = lambda x,t: simplify(gamm1(x,t)*(chi1(x,t)+chi2(x,t)) + chi1(x,t)*chi2(x,t) + 1)
rss = lambda x,t: simplify((-gamm2(x,t)*(chi1(x,t)+chi2(x,t)) - chi1(x,t)*chi2(x,t) + 1)/delta(x,t))
rpp = lambda x,t: simplify((gamm2(x,t)*(chi1(x,t)+chi2(x,t)) - chi1(x,t)*chi2(x,t) + 1)/delta(x,t))
rsp = lambda x,t: simplify(1j*(chi1(x,t) - chi2(x,t))/delta(x,t))
rps = lambda x,t: simplify(-1j*(chi1(x,t) - chi2(x,t))/delta(x,t))
r1 = lambda x,t: np.array([[rss(x,t),rsp(x,t)],[rps(x,t),rpp(x,t)]],dtype='float')
r2 = lambda x,t: r1(x,t)
identity = np.array([[1,0],[0,1]])
g = lambda d,x,t: np.subtract(identity,np.dot(np.dot(r1(x,t),r2(x,t)),identity*e**(-2*x*d)))
funz = lambda d,x,t: h/((2*pi)**2)*np.log(np.linalg.det(g(d,x,t)))*c*sin(t)*x**2
""" ---------------------------------------------- """
def myFunc(x,y):
return funz(10**(-8),x,y)
f = ROOT.TF2("f","myFunc(x,y)",0,10**8,0,2*ROOT.TMath.Pi())
wf1 = ROOT.Math.WrappedMultiTF1(f)
ig = ROOT.Math.AdaptiveIntegratorMultiDim()
ig.SetFunction(wf1)
ig.SetRelTolerance(0.001)
xmin = ROOT.std.vector('double')(2)
xmax = ROOT.std.vector('double')(2)
xmin[0] = 0
xmin[1] = 0
xmax[0] = 10**8
xmax[1] = ROOT.TMath.PiOver2()
print ig.Integral(np.array(xmin),np.array(xmax))
print f.Integral(0,10**8,0,ROOT.TMath.PiOver2())
funz2 = lambda x,y: funz(10**(-8),x,y)
punti = integrate.dblquad(funz2,0,10,0,ROOT.TMath.PiOver2())
print punti
Thanks so much,
Alberto.