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