第四十篇 重复牛顿-柯特斯法则
如果积分的范围很大，一种方法是在整个区间内上拟合一个高阶多项式，然后将这个大区间分解几个等长度的小区间，在每个小区间上使用低阶多项式。上篇数值积分基础-牛顿柯特斯法则中描述的问题都可以使用这种“重复”的形式。
重复矩形法则
如下图所示，将积分区间[A, B]划分为k条，宽度为hi, i = 1,2，…，k。虽然让所有的小区间具有相同的长度是很容易的，但是如果函数变化较快的地方可以让区间更小一些。在一种“自适应”求积的方法中，区间的宽带可以自动修改，这将在后面的文章中描述。
在重复矩形法则中，每个区间的面积近似于一个矩形，其高度由区间下限对应的函数给出。然后把它们加在一起，得到整体的积分值
由上图可以看出，该方法相当于用一系列水平线替换平滑的连续函数f(x)。取得区间越多，越接近函数得实际形状。如果所有的条带宽度都相同h，则公式为
计算实例
使用重复矩阵法则，三条等宽度样条去计算下面积分
得到
因此
精确解为0.7071
可以看出，虽然数值解得精确性较差，但比较上篇数值积分基础-牛顿柯特斯法则单独应用矩阵法则的结果0.5554相比有了明显的改善。而更多的区间可以获得更好的精度。
重复梯形法则
在重复梯形规则中，每个条带的面积近似为一个梯形，如下图所示。
每个梯形在首尾两点上与函数重合。然后把它们加在一起，得到整体的解决方案，因此
由上图可以看出，该方法相当于用一系列线段替换平滑的连续函数f(x)。如果所有的区间都取相同得宽度h，则公式为
计算实例
使用重复梯形法则，3个等分区间来计算下面的积分
得到每段长
因此
精确解为0.7071
可以看出，在这种情况下，通过分成几个区间，数值解得到了很大的改进。因为上篇数值积分基础-牛顿柯特斯法则并未划分区间的单个梯形规则的解只为0.6704。
重复辛普森(simpson)法则
辛普森法则的一整个区间需要三个样本点，因此重复的辛普森法则必须有偶数个小区间，如下图所示。
每对区间必须有相同的宽度，但对与对之间的宽度可以取不同的值。重复的辛普森法则在三个样本点上拟合一条抛物线，假设有k个(偶数)区间，得到以下表达式
如果所有的区间宽度都取h，则规则化简为
计算实例
使用重复辛普森法则，分成四段等宽度，去计算下面积分
每段宽度
精确解为0.7071
在这种情况下，计算结果精确到小数点后4为。在上篇数值积分基础-牛顿柯特斯法则中应用Simpson法则的结果几乎与准确值差不多，为0.7072。
程序如下
其中分为一个主程序，和一个正交牛顿-柯特法则的w和x值的子程序newton_cotes，详情可见上篇数值积分基础-牛顿柯特斯法则的表格
本程序求一个下限为0.25，上限为0.75，重复5次的梯形法则，函数为sinx**2+x的积分值
#重复牛顿-柯特法则
import numpy as np
import math
import B
a=0.25;b=0.75;nsp=2;nr=5
samp=np.zeros((nsp,1))
wt=np.zeros((nsp))
B.newton_cotes(samp,wt)
wr=(b-a)/nr
hr=0.5*wr
area=0
def f61(x):f61=math.sin(x)*math.sin(x)+xreturn f61
for i in range(1,nr+1):cr=a+(i-1)*wr+hrfor j in range(1,nsp+1):area=area+wt[j-1]*hr*f61(cr+samp[j-1,0]*hr)
print('重复牛顿-柯特法则')
print('积分上下限',a,b)
print('牛顿柯特法则的样本点个数',nsp)
print('重复小区间个数',nr)
print('计算结果','{:13.4e}'.format(area))newton_cotes
def newton_cotes(samp,wt):nsp=samp.shape[0]if nsp==1:samp[0,0]=-1.0 wt[0]= 2.0 elif nsp==2:samp[0,0]=-1.0 samp[1,0]= 1.0 wt[0]= 1.0 wt[1]= 1.0 elif nsp==3:samp[0,0]=-1.0 samp[1,0]= 0.0 samp[2,0]= 1.0 wt[0]= 0.33333333333333333333 wt[1]= 1.33333333333333333333 wt[2]= 0.33333333333333333333 elif nsp==4:samp[0,0]=-1.0 samp[1,0]=-0.33333333333333333333 samp[2,0]= 0.33333333333333333333 samp[3,0]= 1.0 wt[0]= 0.25 wt[1]= 0.75 wt[2]= 0.75 wt[3]= 0.25 elif nsp==5:samp[0,0]=-1.0 samp[1,0]=-0.5 samp[2,0]= 0.0 samp[3,0]= 0.5 samp[4,0]= 1.0 wt[0]= 0.15555555555555555556 wt[1]= 0.71111111111111111111 wt[2]= 0.26666666666666666667 wt[3]= 0.71111111111111111111 wt[4]= 0.15555555555555555556 else:print('积分点数量错误')
终端输出结果如下
