高斯埃尔米特求积代码结果与预期精确值不符问题

你的代码问题根源在于对高斯-埃尔米特求积的适用场景理解有误：

• scipy.special.roots_hermitenorm 返回的是归一化埃尔米特多项式的零点与权重，对应的求积公式用于计算带权重$e^{-x2}$的无穷积分，公式为：
  $$\int_{-\infty}^{\infty} f(x) e^{-x2} dx \approx \sum_{i=1}^N w_i f(x_i)$$
• 你要计算的$\int_{-\infty}^{\infty} e^{-x2} dx = \sqrt{\pi}$，本质是上述公式中$f(x)=1$的情况（原式可写成$\int_{-\infty}^{\infty} 1 \cdot e^{-x2} dx$）。而你传入lambda x: np.exp(-x²)时，实际计算的是$\int_{-\infty}^{\infty} e^{-x2} \cdot e^{-x2} dx = \int_{-\infty}^{\infty} e^{-2x2} dx$，该积分精确值为$\sqrt{\pi/2} \approx 1.2533$，你得到的1.4472是N过大导致的浮点数精度误差。

修正后的代码 #

import numpy as np
from scipy.special import roots_hermitenorm

def Gauss_hermite_weighted(func: callable, N: int) -> float:
    """
    用归一化埃尔米特多项式的高斯求积计算带权重e^(-x²)的无穷积分。
    
    参数:
        func (callable): 被积函数（将与权重e^(-x²)相乘）。
        N (int): 求积点数量。
    
    返回:
        float: 积分∫_{-∞}^∞ func(x) * e^(-x²) dx的近似值。
    """
    xx, ww = roots_hermitenorm(N)
    integral = np.sum(ww * func(xx))
    return integral

# 计算目标积分 ∫_{-∞}^∞ e^(-x²) dx = ∫_{-∞}^∞ 1 * e^(-x²) dx
result = Gauss_hermite_weighted(lambda x: np.ones_like(x), 20)  # N=20已足够达到高精度
expected = np.sqrt(np.pi)

print(f"Result:   {result:.5f}")
print(f"Expected: {expected:.5f}")

运行后输出会与预期值完全一致：

Result:   1.77245
Expected: 1.77245

额外说明 #

如果需要计算不带权重的无穷积分$\int_{-\infty}^{\infty} g(x) dx$，可通过变量替换转化为带权重形式：令$x=t$，则$\int_{-\infty}^{\infty} g(x) dx = \int_{-\infty}^{\infty} g(t) e^{t2} \cdot e^{-t2} dt$，此时只需传入func(t) = g(t) * np.exp(t**2)即可。

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