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

Question

我想查看函数 intensity(r) 与 space 的关系图，因此 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

旁白：创建符号表达式时，请保持符号化。不要混合使用真正的浮点数 np.pi 和复杂的浮点数 1j，而是使用 SymPy 的 pi 和 I。

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 似乎不能将贝塞尔函数的乘积与 exp(-x**2) * exp(pi*I*(-x)) 集成。这已经发生在 r 被 1 替换的情况下，因此 r 的符号性质无关紧要。

直接回答你的问题：

Is there any way of performing numerical integration with symbolic variables?

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

Is the only other alternative defining the values of r a priori and finding the integral for each one of them?

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

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

鉴于它是 complex-valued，您还不太清楚您打算如何绘制此函数。也许您是想绘制强度的绝对值？

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

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))

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

Numerical integration vs symbolic integration in python

python

symbolic-math

sympy