将函数与 Python (sympy, quad) 集成，结果是我要绘制的另一个函数答案

【问题标题】：Integrating a function with Python (sympy, quad) where the result is another function I want to plot将函数与 Python (sympy, quad) 集成，结果是我要绘制的另一个函数
【发布时间】：2015-08-06 11:04:00
【问题描述】：

我想使用 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 的积分从零到 x，其中 x 是时间 t 的函数。尝试按原样运行代码会给我一个错误：“ValueError：具有多个元素的数组的真值不明确。使用 a.any() 或 a.all()”。

所以经过短暂的搜索后，我认为也许 sympy 会是要走的路。但是我也失败了。

我想知道是否有人对如何完成这项任务有帮助？

非常感谢，扎克

【问题讨论】：

标签： python scipy sympy integrate quad

【解决方案1】：

您可以对 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() 来完成）。

【讨论】：

【解决方案2】：

我认为您正在寻找的是 lambda 表达式（如果我理解正确的话。请参阅 here 了解更多信息和一些有关 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)

希望我能给你一个体面的回应！

【讨论】：

【解决方案3】：

我认为你得到的错误来自对 scipy.integrate.quad 的错误调用：

第一个参数必须是函数名，然后对该函数的第一个变量执行积分。其他变量的值可以通过 args 关键字传递给函数。
scipy.integrate.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，然后评估表达式 as murison suggested。

【讨论】：