在 Python 中实现复合高斯求积

Question

我想在Python中实现复合高斯求积求积分 ∫₀¹ e^x² dx.使用 Python 的 quad 命令对其进行评估，我得到 ∫₀¹ e^{x ²} dx ≈ 1.46

下面是我在 Python 中实现的尝试。我期望的是随着 n 变大，正交越接近 'real' 积分。然而，当我改变 n 时，结果变得越来越小。我犯了什么错误？

def gauss(f,a,b,n):
    h = float(b-a)/n
    [x,w] = np.polynomial.legendre.leggauss(n)
    result =  0
    for i in range(n):
        result += w[i] * f(x[i])
    result *= h
    return result


    for i in range(1,10):  
        print(gauss(exp,0,1,i))
    
    2.0
    1.3956124250860895
    0.9711551112557443
    0.731135234765899
    0.5850529284514102
    0.4875503086966867
    0.41790049038666144
    0.36566293624426005
    0.32503372130693325

Answer 1

这是基于您的尝试的工作代码：

import numpy as np

def gauss(f,a,b,n):
    half = float(b-a)/2.
    mid = (a+b)/2.
    [x,w] = np.polynomial.legendre.leggauss(n)
    result =  0.
    for i in range(n):
        result += w[i] * f(half*x[i] + mid)
    result *= half
    return result

def fun(x):
    return np.exp(x**2)

for i in range(1,10):  
    print(gauss(fun,0,1,i))

问题是积分必须重新缩放到您使用的勒让德多项式的域。此外，您没有正确使用正交规则（它没有使用代码中的步长 h。）正确的规则是：

∫_a^bf(x) dx ≈ (b-a)/2 ∑_i=1ⁿw_i f((b-a)/2ξ_i + (a+b)/2).

上面的代码输出：

1.2840254166877414
1.4541678892391303
1.462409711477322
1.462646815656646
1.4626516680186825
1.4626517449041974
1.4626517458962964
1.4626517459070796
1.462651745907181

在 Python 中实现复合高斯求积

Implementing composite Gauss quadrature in Python

python

numerical-methods