python中的数值积分与符号积分

Question

我想看一个函数intensity(r)与空间的图，所以r（径向对称）。

但是，我的强度来自intensity(r) = integrate(integrand(r), (x,0,5))，其中integrand = exp(-x**2) * exp(np.pi*1j*(-x)) * besselj(0, r*x) * x。

上面的所有语法都使用了sympy包，所以我首先定义了 x,y = symbols('x r') .

---

我使用符号变量是因为它似乎使视觉上更容易，r直到最后，当我绘制它并为其分配一个数值时，它才作为变量离开。

然而，用符号变量做那个可怕的积分似乎需要很长时间。

• 有没有办法用符号变量进行数值积分？

• 定义r先验值并为每个值找到积分的唯一另一种选择是什么？

Answer 1

旁白：当你创建一个象征性的表达时，保持它的象征性。不要混合在一个真正的 floatnp.pi和一个复杂的 float 中1j，使用 SymPy's piandI代替。

from sympy import exp, pi, I, besselj, symbols
x, r = symbols('x r')
integrand = exp(-x**2) * exp(pi*I*(-x)) * besselj(0, r*x) * x

但是，是的，SymPy 似乎不能将 Bessel 函数的乘积与exp(-x**2) * exp(pi*I*(-x)). 这已经发生在 r 被 1 替换的情况下，所以 r 的符号性质是无关紧要的。

直接回答你的问题：

有没有办法用符号变量进行数值积分？

没有，就像没有干水一样。这在术语上是矛盾的。

唯一的另一种选择是定义先验值并为每个值找到积分吗？

是的。它可以通过 SymPy（将调用 mpmath）完成：

>>> intensity = lambda r_: Integral(integrand.subs(r, r_), (x, 0, 5)).evalf()
>>> intensity(3)
0.0783849036516177 - 0.125648626220306*I

考虑到它是复数值，您计划如何绘制此函数还不是很清楚。也许您打算绘制强度的绝对值？

无论如何，与 SymPy/mpmath（纯 Python）集成对于绘图来说太慢了。您最好使用 SciPyquad进行集成。它不处理复杂的被积函数，所以我分别对实部和复数部分进行积分。

from scipy.integrate import quad
from scipy.special import jn
integrand = lambda x, r: np.exp(-x**2) * np.exp(np.pi*1j*(-x)) * jn(0, r*x) * x
intensity = lambda r: np.sqrt(quad(lambda x: np.real(integrand(x, r)), 0, 5)[0]**2 + quad(lambda x: np.imag(integrand(x, r)), 0, 5)[0]**2)

现在intensity(3)的评估速度比以前的版本快得多。我们可以绘制它：

import matplotlib.pyplot as plt
t = np.linspace(0, 3)
plt.plot(t, np.vectorize(intensity)(t))