我想使用python对一个函数进行积分,其中输出是一个新函数而不是一个数值。例如,我有一个方程(来自Arnett 1982——超新星的解析描述):
def A(z,tm,tni):
y=tm/(2*tni)
tm=8.8 # diffusion parameter
tni=8.77 # efolding time of Ni56
return 2*z*np.exp((-2*z*y)+(z**2))
我想求出A的积分,然后画出结果。首先,我天真地尝试了scipy.quad:
def Arnett(t,z,tm,tni,tco,Mni,Eni,Eco):
x=t/tm
Eni=3.90e+10 # Heating from Ni56 decay
Eco=6.78e+09 # Heating from Co56 decay
tni=8.77 # efolding time of Ni56
tco=111.3 # efolding time of Co56
tm=8.8 # diffusion parameter
f=integrate.quad(A(z,tm,tni),0,x) #integral of A
h=integrate.quad(B(z,tm,tni,tco),0,x) #integral of B
g=np.exp((-(x/tm)**2))
return Mni*g*((Eni-Eco)*f+Eco*h)
其中B也是一个预定义函数(这里没有给出)。A和B都是z的函数,但是最后的方程是时间t的函数(我相信正是在这里我导致了我的代码失败)
A和B的积分从0到x,其中x是时间t的函数。尝试按照它的形式运行代码会给我一个错误:"ValueError:包含多个元素的数组的真值是模糊的。使用a.a any()或a.a all()".
所以经过简短的搜索,我想也许sympy会是一种方法。然而,我在这方面也失败了。
我想知道是否有人有一个有用的建议如何完成这项任务,请问?很多谢谢,扎克
可以解析地积分A。假设我没有因为起床太晚而错过一些愚蠢的事情,下面的建议对你有帮助吗?
import sympy as sy
sys.displayhook = sy.pprint
A, y, z, tm, t, tni = sy.symbols('A, y, z, tm, t, tni')
A = 2*z*sy.exp(-2*z*y + z**2)
expr = sy.integrate(A, (z,0,t)) # patience - this takes a while
expr
# check:
(sy.diff(expr,t).simplify() - A.replace(z,t)).simplify()
# thus, the result:
expr.replace(y,tm/(2*tni)).replace(t,t/tm)
最后一行生成解析形式的A函数的积分,尽管它确实需要计算假想的误差函数(可以使用scipy.special.erfi())。
我认为你正在寻找的是lambda表达式(如果我正确理解你所说的…有关lambda函数的更多信息和一些示例,请参阅此处。
他们允许你做的是在A中定义一个匿名函数并返回它,这样你就得到了B函数,应该是这样的:
def A(parameters):
return lambda x: x * parameters # for simplicity i applied a multiplication
# but you can apply anything you want to x
B = A(args)
x = B(2)
希望我能给你一个像样的答复!
我认为您得到的错误来自于对scipy. integration .quad的错误调用:
- 第一个参数需要只是函数名,然后对该函数的第一个变量执行积分。其他变量的值可以通过args关键字传递给函数。
- scipy. integrated .quad的输出不仅包含积分值,还包含误差估计。所以返回一个包含2个值的元组!
最后,下面的函数应该可以工作:
def Arnett(t, z, Mni, tm=8.8, tni=8.77, tco=111.3, Eni=3.90e+10,
Eco=6.78e+09):
x=t/tm
f,err=integrate.quad(A,0,x,args=(tm,tni)) #integral of A
h,err=integrate.quad(B,0,x,args=(tm,tni,tco)) #integral of B
g=np.exp((-(x/tm)**2))
return Mni*g*((Eni-Eco)*f+Eco*h)
但更好的解决方案可能是对A和B进行分析积分,然后像murison建议的那样对表达式求值。