当前位置：首页 > article >正文

Python数值计算（32）——simpson 1/3积分公式

article 2024/10/25 7:53:12

1. 背景知识

前面我们通过用矩形和梯形的数值算法，近似实现了数值积分，那么，和之前插值类似，是否可以使用多项式来拟合曲线，然后将该多项式作为被积函数求积分呢？当然是可行的，如果以最简单的二次多项式为例，如果给出三个点及对应的函数值，一般就可以构造出一个二次多项式，然后以此进行积分，这就是Simpson（辛普森）积分公式的思想。

通过已知的三对数 $(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))$ ，采用Lagrange插值法，可以得到如下公式：

$\int_{a}^{b}f(x)dx=\int_{a}^{b}\sum_{i=0}^{n}\prod_{j=0,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}f(x_i)dx$

简化一下形式可以写成：

$\int_{a}^{b}f(x)dx=\int_{a}^{b}(px^2+qx+r)dx=(px^3/3+qx^2/2+rx)|_{a}^{b}$

其中 $p,q,r$ 是待定系数，如果我们取均匀的点，例如 $x_0=a,x_1=(a+b)/2,x_2=b$ ，则

可得上式为：

$(\frac{px^3}{3}+\frac{qx^2}{2}+rx)|_{a}^{b}=\frac{h}{3}(f(x_0)+4f(x_1)+f(x_2))$

其中：

$h=\frac{b-a}{2}$

由于前面的系数是1/3，因此得名辛普森1/3积分公式。

2. Simpson 1/3积分

掌握了上述原理以后，实现Simpson 1/3积分公式将十分容易：

# Simpson's 1/3 Rule for Numerical Integration
import numpy as np
def simpson13(f,a,b):
    """
    Simpson's 1/3 Rule for Numerical Integration
    """
    h=(b-a)/2
    x=np.array([a,a+h,b])
    y=np.vectorize(f)(x)
    return h/3*(y[0]+4*y[1]+y[2])

和前面的算法类似，单个区间的误差还是比较大的：

3. 复合Simpson 1/3积分公式

同样的，单个区间的误差较大，我们可以将一个大的区间划分为多个小区间，然后分别在 $[x_0,x_2],[x_2,x_4],...,[x_{2i},x_{2i+2}],...,[x_{n-2},x_n]$ 的小区间上进行Simpson 1/3积分公式，这也就是复合Simpson 1/3积分公式。

算法实现为：

import numpy as np

def comp_simpson13(f,a,b,n):
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    y=np.vectorize(f)(x)    
    return h/3*(y[0]+4*sum(y[1:n:2])+2*sum(y[2:n:2])+y[n])

多个区间的效果如下：