Python中数值积分与符号积分对比及径向强度函数绘图问题

我来帮你捋捋这个问题哈——你想用Sympy定义符号变量来保留r，最后绘图时赋值，但遇到了麻烦对吧？核心问题在于你这个积分没有解析解，Sympy没法直接给出一个关于r的闭合表达式，所以直接用符号积分后再代入数值会要么卡很久，要么报错。下面给你两个实用的解决方案：

方案一：用Scipy数值积分（推荐，效率更高） #

既然最终是要绘图出数值结果，直接用Scipy的数值积分工具会更高效，毕竟Sympy擅长符号推导，数值计算还是Scipy更拿手。

import numpy as np
from scipy.integrate import quad
from scipy.special import j0  # 第一类零阶贝塞尔函数的数值实现
import matplotlib.pyplot as plt

def calculate_intensity(r):
    # 定义被积函数，这里x是积分变量，r是输入参数
    def integrand(x):
        exp_term1 = np.exp(-x**2)
        exp_term2 = np.exp(np.pi * 1j * (-x))  # 复指数项
        bessel_term = j0(r * x)
        return exp_term1 * exp_term2 * bessel_term * x
    
    # 复值积分需要分开计算实部和虚部，Scipy的quad默认处理实函数
    real_part, _ = quad(lambda x: np.real(integrand(x)), 0, 5)
    imag_part, _ = quad(lambda x: np.imag(integrand(x)), 0, 5)
    return real_part + 1j * imag_part

# 生成r的取值范围，这里取0到10，100个点
r_values = np.linspace(0, 10, 100)
# 计算每个r对应的强度模值（如果需要实部/虚部直接改就行）
intensity_magnitudes = np.array([np.abs(calculate_intensity(r)) for r in r_values])

# 绘图
plt.figure(figsize=(8, 5))
plt.plot(r_values, intensity_magnitudes, color='#1f77b4', linewidth=2)
plt.xlabel('r')
plt.ylabel('|intensity(r)|')
plt.title('Radially Symmetric Intensity vs. Spatial Variable r')
plt.grid(alpha=0.3)
plt.show()

方案二：坚持用Sympy的符号+数值混合方案 #

如果你确实需要先保留符号形式做一些推导，之后再转数值绘图，可以用lambdify把Sympy的积分表达式转成可计算的数值函数：

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# 定义符号变量
x, r = sp.symbols('x r')
# 定义被积函数（注意用Sympy的内置函数，别混numpy的）
integrand = sp.exp(-x**2) * sp.exp(sp.I * sp.pi * (-x)) * sp.besselj(0, r*x) * x

# 把积分表达式转成数值函数，指定用numpy和scipy做计算后端
intensity_func = sp.lambdify(
    r, 
    sp.integrate(integrand, (x, 0, 5), conds='none'),
    modules=['numpy', 'scipy']
)

# 生成r的数值（这里点数别太多，Sympy转的函数速度会慢一点）
r_values = np.linspace(0, 10, 50)
# 计算强度模值
intensity_magnitudes = np.abs(intensity_func(r_values))

# 绘图
plt.figure(figsize=(8, 5))
plt.plot(r_values, intensity_magnitudes, color='#ff7f0e', linewidth=2)
plt.xlabel('r')
plt.ylabel('|intensity(r)|')
plt.title('Intensity vs. r (Sympy-based)')
plt.grid(alpha=0.3)
plt.show()

为啥原来的方法会出问题？ #

你定义的被积函数结合了高斯函数、复指数和贝塞尔函数，这类积分没有解析形式的解，Sympy的sp.integrate只会返回一个未计算的Integral对象，没法直接代入r的数值得到结果。必须通过数值积分的方式来计算每个r对应的强度值。

内容的提问来源于stack exchange，提问作者SuperCiocia