我想将包含符号变量的表达式转换为数字表达式,以便该表达式随后可以在积分方法"quad"中使用。
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(ktsq,qbarsq,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))
代码运行到这里并返回正确的函数fA
。fA
是一个符号值表达式,我想传递给 quad 以执行积分(超过ktsq
)
import scipy.integrate.quadrature as sciquad
def integrated_f(qbarsq,NN,a,b,w):
return sciquad(fA,1,(w**2-MJpsi**2)/4, args=(qbarsq, NN, a, b, w))
我的理解是,这失败了,因为quad
的第一个参数,即函数,是符号类型,而不是四边形所需的数字(=浮点)。如何使函数数字化,从而允许我执行集成?我已经尝试了.subs
和lambdify
功能,但无法使其工作。前者似乎仅在提供数字时才有效(例如,我不想这样做NN=0.1
设置),并且我尝试了以下内容进行lambdify
def test(ktsq):
return fA(ktsq,qbarsq,NN,a,b,w)
f = lambdify(((qbarsq,NN,a,b,w),), test(ktsq))
#print(f(1,2,3,4,5))
但是当我取消注释打印以检查是否一切正常时,这给出了有关位置参数数量的错误。
TypeError: <lambda>() takes 1 positional argument but 5 were given
是的,你应该使用 lambdify。lambdify
的第一个参数是符号元组,而不是代码中的元组元组。第二个参数是 SymPy 表达式。例:
from sympy import *
a, b, c, d, e = symbols('a b c d e')
expr = a*b + 2*c + d/e
f = lambdify((a, b, c, d, e), expr)
print(f(1, 2, 3, 4, 5)) # prints 8.8
在您的情况下,这看起来像
expr = fA(qbarsq, ktsq, NN, a, b, w)
f = lambdify((qbarsq, ktsq, NN, a, b, w), expr, "mpmath")
这里选择 mpmath 作为后端,因为它可以计算表达式包含的 Gamma 函数。否则,使用"numpy"后端选项可能会更快。查看更多关于lambdify的信息。
print(f(1, 2, 3, 4, 5, 6)) # (-4757.21371513605 + 58978.7828908493j)
它如何与quad
一起工作...取决于你得到的是实数还是复数。当它们是真实的时,您可以与scipy.integrate.quad
集成:
from scipy.integrate import quad
quad(f, 3, 4, args=(2, 3, 4, 5, 6))[0]
返回30049812.82526324
。
如果它们很复杂,SciPy 的quad
会对它不理解mpc
类型不满意。但是mpmath
有自己的quad
,所以改用它:
import mpmath as mp
mp.quad(lambda x: f(2, x, 3, 4, 5, 6), [1, 3])
返回mpc(real='7170810.3848631922', imag='-192389955826656.31')
。这里 [1, 3] 是积分的区间。