I have defined a series of functions below with the end function fA
being inserted for a numerical integration. The integration is with respect with one variable so the other arguments are passed as numeric so that the integration method (quad
) may proceed
import numpy
import math as m
import scipy
import sympy
#define constants
gammaee = 5.55e-6
MJpsi = 3.096916
alphaem = 1/137
lambdasq = 0.09
Ca = 3
qOsq = 2
def qbarsq(qsq):
return (qsq+MJpsi**2)/4
def xx(qbarsq, w):
return 4*qbarsq/(4*qbarsq-MJpsi**2+w**2)
from sympy import *
x,NN,a,b,ktsq,qbarsq,w = symbols('x NN a b ktsq qbarsq w')
def xg(a,b,NN,ktsq,x):
return NN*(x**(-a))*(ktsq**b)*exp(sqrt((16*Ca/9)*log(1/x)*log((log(ktsq/lambdasq))/(log(qOsq/lambdasq)))))
#prints symbolic derivative of xg
def func(NN,a,b,x,ktsq):
return (-x*diff(log(xg(a,b,NN,ktsq,x)),x))
#print(func(NN,a,b,x,ktsq))
#prints symbolic expression for Rg
def Rg(NN,a,b,ktsq,x):
return 2**(2*func(NN,a,b,x,ktsq)+3)/sqrt(m.pi)*gamma(func(NN,a,b,x,ktsq)+5/2)/gamma(func(NN,a,b,x,ktsq)+4)
#print(Rg(NN,a,b,ktsq,x))
#prints symbolic expression for Fktsq
def FktsqDeriv(NN,a,b,x,ktsq):
return diff(Rg(NN,a,b,ktsq,x)*xg(a,b,NN,ktsq,x),ktsq)
#print(FktsqDeriv(NN,a,b,x,ktsq))
def Fktsq1(qbarsq,ktsq,NN,a,b,w):
return FktsqDeriv(NN,a,b,x,ktsq).subs(x,4*qbarsq/(4*qbarsq-MJpsi**2+w**2))
print(Fktsq1(qbarsq,ktsq,NN,a,b,w))
# symbolic expression for fA
def fA(qbarsq,ktsq,NN,a,b,w):
return Fktsq1(qbarsq,ktsq,NN,a,b,w)*1/(qbarsq)*1/(qbarsq+ktsq)
#print(fA(qbarsq,ktsq,NN,a,b,w))
I now want to integrate this last function over ktsq
as follows,
import scipy.integrate.quadrature as sciquad
def integrated_f(NN,a,b,w,qbarsq):
return sciquad(fA,1,(w**2-MJpsi**2)/4, args=(NN, a, b, w, qbarsq))
a=0.1
NN=0.5
b=-0.2
w=89
qbarsq=5
result = integrated_f(NN,a,b,w,qbarsq)
print(result)
The problem is where I try to get a number out from this integration by specifying numerical values for each of the other parameters. The error is
ValueError:
Can't calculate 1st derivative wrt 989.426138911118.
My only interpretation of this is that the method cannot cope with the complexity of the function (even though I think it is relatively simple in structure) because I did not define any more derivatives and certainly not with respect to this value. Is there an easy solution? Actually I wish to use the function integrated_f
to use in an optimisation problem for best fit parameters a,b,NN
. Would something like
scipy.optimize.minimize(integrated_f, x0, method='Nelder-Mead', options={'max\
iter': 1000})
be ok for multivariable functions where x0
is an array of initial guesses. Thanks!