高斯求积与埃尔米特多项式的问题

Question

我写了一段代码，但结果和我预期的差得很远，我本来希望得到一个更接近正确值的结果。

import numpy as np
from scipy.special import roots_hermitenorm

def Gauss_hermite(func: callable, N: int) -> float:
    """
    Numerically computes the integral of a function using Gauss quadrature
    with Hermite polynomials.

    Parameters:
        func (callable): The function to be integrated.
        N (int): The number of intervals for integration.

    Returns:
        float: Approximate value of the integral of 'func' over the interval
            [a, b].
    """

    # Determine zeros and weights using scipy
    xx, ww = roots_hermitenorm(N)
    
    # Compute the integral
    integral = np.sum(ww * func(xx))
    
    return integral

# Example usage
result = Gauss_hermite(lambda x: np.exp(-x**2), 2046)
expected = np.sqrt(np.pi)

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

结果是：

实际结果：1.44720
预期结果：1.77245

Answer 1

根据维基百科上的高斯-赫尔米特积分法，从概念上讲，你需要的东西大概是这样的：

integral = np.sum(ww * func(xx)/np.exp(-xx**2/2))

这个积分公式是用来计算一个被np.exp(-xx**2/2)加权的函数（因为SciPy文档提到这些多项式是正交的，权重函数是np.exp(-x**2/2)，而不是np.exp(-x**2)），所以你需要把这个加权去掉。

对于较低阶的多项式（比如64），这样做能得到合理的结果，但如果阶数像2048那么高，就会遇到数值问题。所以，实际上，与其改变权重，不如通过数学方式把你的被积函数除以np.exp(-x**2/2)来改变它：

result = Gauss_hermite(lambda x: np.exp(-x**2/2), 2046)

如果你有另一个被积函数，无法那么简单地去掉加权，那么可以用其他方法来解决数值问题，或者使用更合适的积分规则，但那是另一个问题。

---

改变权重：

import numpy as np
from scipy.special import roots_hermitenorm

def Gauss_hermite(func: callable, N: int) -> float:
    xx, ww = roots_hermitenorm(N)
    return np.sum(ww * func(xx)/np.exp(-xx**2/2))

res = Gauss_hermite(lambda x: np.exp(-x**2), 64)
print(res)  # 1.7724538509055154
np.testing.assert_allclose(res, np.sqrt(np.pi))  # passes

改变被积函数：

def Gauss_hermite(func: callable, N: int) -> float:
    xx, ww = roots_hermitenorm(N)
    return np.sum(ww * func(xx))

res = Gauss_hermite(lambda x: np.exp(-x**2/2), 2046)
print(res)  # 1.7724538509055427
np.testing.assert_equal(res, np.sqrt(np.pi))  # passes