1. 背景知识
前面我们通过用矩形和梯形的数值算法,近似实现了数值积分,那么,和之前插值类似,是否可以使用多项式来拟合曲线,然后将该多项式作为被积函数求积分呢?当然是可行的,如果以最简单的二次多项式为例,如果给出三个点及对应的函数值,一般就可以构造出一个二次多项式,然后以此进行积分,这就是Simpson(辛普森)积分公式的思想。
通过已知的三对数,采用Lagrange插值法,可以得到如下公式:
简化一下形式可以写成:
其中是待定系数,如果我们取均匀的点,例如,则
可得上式为:
其中:
由于前面的系数是1/3,因此得名辛普森1/3积分公式。
2. Simpson 1/3积分
掌握了上述原理以后,实现Simpson 1/3积分公式将十分容易:
python"># 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)/2x=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积分公式
同样的,单个区间的误差较大,我们可以将一个大的区间划分为多个小区间,然后分别在的小区间上进行Simpson 1/3积分公式,这也就是复合Simpson 1/3积分公式。
算法实现为:
python">import numpy as npdef comp_simpson13(f,a,b,n):h=(b-a)/nx=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])
多个区间的效果如下: