从Newton-Cotes的截断误差公式可以看出，当积分区间$[a,b]$较大时，低阶的Newton-Cotes求积公式截断误差都比较大。由于高阶Newton-Cotes求积公式是数值不稳定的，因此通过不断增加阶数来提高求积公式的精度是不可行的。但是，如果将积分区间$[a,b]$分成几个小区间（任意的），在每个小区间上应用Newton-Cotes求积公式，其截断误差必然会减小，然后再把每个小区间上的积分值累加起来，这样却能大大提高整个积分的精度。这种方法称为复化求积方法。
常用的复化求积方法采用等分区间的做法，具体如下：
将区间$[a,b]$划分为n等分，步长为$H=\frac{(b-a)}{n}$，分点为$x_k=a+kH,k=0,1,2,\cdots,n$。先用低阶Newton-Cotes求积公式求得每个子区间$[x_k,x_{k+1}]$上的积分值$I_k$，然后将它们累加起来求和，用$\sum_{k=0}^{n-1}I_k$作为所求积分$I=\int_a^bf(x)dx$的近似值。
1. 复化梯形公式
在区间$[a,b]$上采用复化求积方法，具体使用梯形求积公式进行计算，就得到复化梯形求积公式。用$T_k$表示$f(x)$在子区间$[x_k,x_{k+1}]$上的积分值，$T_n$表示$f(x)$在区间$[a,b]$上的积分值，有：

$T_k=\frac{1}{2}H[f(x_k)+f(x_{k+1})]$

其中

$H=(b-a)/n, \quad x_k=a+kH \quad(k=0,1,2,\cdots,n)$

故

$T_n=\sum_{k=0}^{n-1}T_k=\frac{1}{2}H\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]$

即：

$T_n=\frac{1}{2}H[f(x_{k=0})+f(x_{k=n})+2\sum_{k=1}^{n-1}f(x_k)]$

或

$T_n=\frac{b-a}{2n}[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_k)]$

截断误差用$R_T$表示，由于$f(x)$在区间$[a,b]$上有连续的二阶导数，故有：

$R_T=\sum_{k=0}^{n-1}-\frac{1}{12}H^3f^{(2)}(\eta_k)=-\frac{1}{12}H^3\sum_{k=0}^{n-1}f^{(2)}(\eta_k) \quad \eta_k\in[x_k,x_{k+1}]$

即

$R_T=-\frac{1}{12}H^3·n·f^{(2)}(\eta)=-\frac{(b-a)}{12}H^2f^n(\eta) \quad \eta \in [a,b]$
2. 复化Simpson公式
在区间$[a,b]$上采用复化求积方法，具体使用Simpson求积公式进行计算，就得到复化Simpson公式。用$S_k$表示$f(x)$在子区间$[x_k,x_{k+1}]$上的积分值，$S_n$表示$f(x)$在区间$[a,b]$上的积分值，有：

$S_k=\frac{1}{6}H[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$

式中，$x_{k+\frac{1}{2}}$为子区间$[x_k,x_{k+\frac{1}{2}}]$的中点，$H=(b-a)/n$。

$S_n=\sum_{k=0}^{n-1}=\frac{1}{6}H\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$

故

$S_n=\sum_{k=0}^{n-1}=\frac{1}{6}H[f(a)+4\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum_{k=1}^{n-1}f(x_{k})+f(b)]$

截断误差用$R_s$表示，由于$f(x)$在区间$[a,b]$上有连续的四阶导数，故有：

$R_S=\sum_{k=0}^{n-1}-\frac{1}{2880}H^5f^{(4)}(\eta_k)=-\frac{1}{2880}H^5\sum_{k=0}^{n-1}f^{(4)}(\eta_k) \quad \eta_k \in [x_k,x_{k+1}]$
3. 复化Cotes公式
在区间$[a,b]$上采用复化求积方法，具体使用Cotes求积公式进行计算，就得到复化Cotes公式为：

$C_n=\frac{1}{90}H[7f(a)+32\sum_{k=0}^{n-1}f(x_{k+\frac{1}{4}})+12\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+32\sum_{k=0}^{n-1}f(x_{k+\frac{3}{4}})+14\sum_{k=1}^{n-1}f(x_k)+7f(b)]$

截断误差为：

$R_c=-\frac{2(b-a)}{945}(\frac{H}{4})^6f^{(6)}(\eta) \quad \eta\in[a,b]$

从复化求积的余项公式中可以看出，复化梯形公式、复化Simpson公式、复化Cotes公式的余项和步长的关系为$R_T=O(h^2),R_s=O(h^4),R_c=O(h^6)$。因此，当$H\to 0$或$n\to \infty$时，$T_n,S_n,C_n\to I$。

对于函数f(x) = sin(x) ，试用复化梯形公式、复化Simpson公式及复化Cotes公式计算积分，并比较其误差。

首先，在MATLAB输入以下程序定义函数：

function y=f(x) 
y=sin(x);


1.复化梯形求积公式的代码

function Tn=Tn(a,b,n)
    format long
    h=(b-a)/n;
    sum=0;
    for k=1:n-1
        sum=sum+f(a+k.*h);
    end
    Tn=(f(a)+2*sum+f(b))*h/2;
  end



2.复化Simpson公式的代码

function Sn = Sn(a,b,n)
    format long
    h = (b-a)/n;
    sum1 = 0;
    sum2 = 0;
    for i = 0:n-1
        sum1 = sum1 + f(a+(i+1/2).*h);
    end
    for j = 1:n-1
        sum2 = sum2 + f(a+j.*h);
    end
    Sn = h/6*(f(a)+4*sum1+2*sum2+f(b));



3.复化Cotes公式的代码

function Cn = Cn(a,b,n)
    format long
    h = (b-a)/n;
    sum1 = 0;
    sum2 = 0;
    for i = 0:n-1
        sum1 = sum1 + 32*f(a+(i+1/4).*h)+12*f(a+(i+1/2).*h)+32*f(a+(i+3/4).*h);
    end
    for j = 1:n-1
        sum2 = sum2 + 14*f(a+j.*h);
    end
Cn = h/90*(7*f(a)+sum1+sum2+7*f(b));



4.不同区间和分割次数下的精确值和近似值的比较代码

y=[-0.150920964949986,0.106044306328975,-1.260000791649389,1.122733727670800,0.185427920752058,0.955618503172607];
Tn=[-0.144707438957434 ,0.103357516942886,-0.809037206720437,1.074587530720150,0.176653381231435 ,0.874739367319743];
Sn=[-0.150933735388946,0.106047691917755 ,-1.267931587831384,1.122836754412495,0.185448626084829 ,0.956086587141625];
Cn=[ -0.150920958200157,0.106044290463326,-1.259953005777527,1.122733636186552,0.185427890841067,0.955744078179222];
plot(y,'-k*');
hold on;
plot(Tn,'-g.');
hold on;
plot(Sn,'-bo');
hold on;
plot(Cn,'-r+');
legend('精确值','Tn','Sn','Cn');



5.图像


6.结果分析
1 误差分析

在上面的数值计算中，复化Cotes公式计算所得的误差最小，其次是复化Simpson公式，误差最大的是复化梯形公式。


2 精确值比较

它们的精度有很大的差别：与精确值比较，复化梯形公式的结果只有一位有效数字，复化Simpson公式的结果有四位有效数字，而复化Cotes公式的结果却有七位有效数字。

由此可知，复化Cotes公式的代数精度最高，复化Simpson其次，复化梯形公式的代数精度最低。


3 收敛性分析



4 稳定性分析



可知三个求积公式的求积系数都为正，所以复化梯形公式、复化Simpson公式和复化Cotes公式均是稳定的。

由误差余项公式可知区间过大，误差亦大；为避免可选取适当多的节点，即选取相对高阶的Newton-cotes公式，但由稳定性分析又知：当阶数过大时，会出现不稳定的Runge现象。
复化求积法：将积分区间进行适当分段，在各分段子区间上采用低阶的Newton-Cotes求积公式，对各个小区间上的积分值进行一个近似，最后再累加起来。
例如：复化Simpson公式的推导：
复化Simpson公式误差分析：
其中有个加一项减一项 
区间逐次二分法
由复化求积公式的截断误差可知，加密节点可以提高求积公式的精度，但困难在于：使用公式之前需给出合适的步长，h过大，满足不了精度；h过小，计算量过大，因而实用的方法是采用区间逐次二分法，反复利用求积公式计算，直至二分前后两次积分值的差满足精度为止。
比如：对区间进行n等分，每个区间上先采用梯形公式，即复化梯形公式，若不能满足精度，则将每个小区间二等分，再分别采用梯形公式，不过端点处的值不用再算了，新算的只有新小区间上的二等分点处值，这样便可使计算量节约一半。
龙贝格算法（自动调整等分数）

实
验
目
的
或
要
求
1
、利用复化梯形公式、复化
simpson
公式计算积分
2
、比较计算误差与实际误差
实
验
原
理
(
算
法
流
程
图
或
者
含
注
释
的
源
代
码
)
取
n=2,3,
…
,10
分别利用复化梯形公式、复化
simpson
公式计算积分
1
2
0
I
x
dx


，
并与真值进行比较，
并画出计算误差与实际误差之间的曲线。
利用复化梯形公式的程序代码如下：
function f=fx(x)
f=x.^2;
％首先建立被积函数，以便于计算真实值。
a=0;
％积分下线
b=1;
％积分上线
T=[];
％用来装不同
n
值所计算出的结果
for n=2:10;
h=(b-a)/n;
％步长
x=zeros(1,n+1);
％给节点定初值
for i=1:n+1
x(i)=a+(i-1)*h;
％给节点赋值
end
y=x.^2;
％给相应节点处的函数值赋值
t=0;
for i=1:n
t=t+h/2*(y(i)+y(i+1));
％利用复化梯形公式求值
end
T=[T,t];
％把不同
n
值所计算出的结果装入
T
中
end
R=ones(1,9)*(-(b-a)/12*h.^ 2*2);
％积分余项(计算误差)
true=quad(@fx,0,1);
％积分的真实值
A=T-true;
％计算的值与真实值之差(实际误差)
x=linspace(0,1,9);
plot(x,A,'r',x,R,'*')
％将计算误差与实际误差用图像画出来
注：由于被积函数是
x.^2
，它的二阶倒数为
2
，所以积分余项为：
(-(b-a)/12*h.^ 2*2)