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

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)

print("Result of integration:", result)

输出：1.4472025091165737 精确值：sqrt(pi) ~ 1,77245

Answer 1

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

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

求积公式用于评估由

np.exp(-xx**2/2)

加权的被积函数，因此您需要撤消该加权。

这为低阶多项式提供了合理的结果，但您会遇到像

2048

这样的阶数的数值问题。