厄米多项式的高斯求积问题

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

基于维基百科上的Gauss-Hermite Quadrature，看起来（概念上）你想要类似的东西：

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