第6章：数值积分与数值微分
6.1 求积公式由定积分的定义可知，连续函数f(x) 在区间[a,b]上的定积分近似值可以表示为[a,b]内的一些点x0,x1,…,xn处的函数值f(x0),f(x1),…,f(xn)的加权和或线性组合，即
（1）
其中w0,w1,…,wn仅与x0,x1,…,xn有关而与被积函数f(x)无关。我们把这样的公式称为求积公式，也称为机械求积公式。
1．术语和记号为了计算f(x) 在区间[a,b]上的定积分近似值，我们通常的做法是，把积分区间[a,b]划分为n等分，记h=(b-a)/n,x0=a,xk=a+kh,k=0,1,2,…,n,称x0,x1,…,xn为[a,b]的一个等份分划。
假如x0,x1,…,xn为[a,b]的一个等份分划那么求积公式(1)中的w0,w1,…,wn的选取仅仅只与n有关，从而可以简化对求积公式的研究。
2.求积公式的性质微积分学中我们曾研究过，定积分保持函数的线性关系不变，它的含义是，若f(x),g(x)都是[a,b]上的可积函数，则对任意实数u,v,我们有u·f(x)+v·g(x)也是[a,b]上的可积函数，而且
不难验证，求积公式也保持函数的线性关系不变，即
3.几种常见的求积公式在后面的讨论中，我们将经常用到下面一些非常简单的求积公式，他们是中点公式、梯形公式和辛卜生公式。
4．截断误差在求积公式中，我们使用的是近似等号，这是因为，对于一般的被积函数来说，利用这些公式计算所得的结果除了舍入误差外，还有截断误差。
有时为了进行误差分析，我们可以把求积公式式写成
（1’）
其中R[f]表示的就是截断误差。
5.代数精度的概念定义：一个求积公式
如果对所有的次数不超过m的多项式严格相等，而对某些m+1次多项式不相等，则称该公式具有代数精度m，或该公式的代数精度为m。
利用求积公式的线性性，我们不难证明下面的结论。
定理：如果求积公式对1，x,…,xm严格相等，而对xm+1不相等，则该公式的代数精度为m。
作为课外练习，鼓励大家给出完整证明。
6.基本结论我们不难利用上面的定理所给出的方法证明辛卜生公式的代数精度是3，而中点公式和梯形公式的代数精度是1。
现在我们可以对这三个公式作一个简单的评价：
中点公式和梯形公式的代数精度虽然都是1，但中点公式只计算一个点的函数值，而梯形公式要计算两个点处的函数值，所以中点公式优于梯形公式。
与梯形公式相比，辛卜生公式只多计算一个点的函数值，但代数精度却增加到3，显然辛卜生公式更为优越。
6．2 牛顿-柯特斯求积公式利用插值多项式近似替代被积函数设f(x)为被积函数，[a,b]为积分区间，x0,x1,…,xn为[a,b]内的n+1个互异的点，记Ln(x)为相应的拉格朗日插值多项式，那么我们有
2.利用插值多项式导出求积公式
在上面给出的公式中，由于诸lk(x)都是多项式函数，所以诸wk都可以精确地计算出来。从而我们可以得到一般性的求积公式。
3.讨论：
回顾上一章关于多项式插值的结论，由于任意次数不超过n的多项式与它的任意n+1个基点的插值多项式恒等，再由求积公式的代数精度的定义，我们立即得到：由n个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少式n。
为了同时说明求即公式得代数精度，我们记
4.牛顿-柯特斯求积公式牛顿-柯特斯求积公式就是利用等距基点的拉格朗日插值多项式导出的求积公式。
将积分区间[a,b]划分为n等分，记h=(b-a)/n，取
x0=a,xk=a+kh.k=0,1,…,n
我们可以得到
在牛顿-柯特斯公式中，我们称c(n,k)为牛顿-柯特斯系数，一般可通过查表得到。崔国华教材p60列出了直到n=8的所有牛顿-柯特斯系数，应该说，实用意义不大。崔国华教材p59-60给出了牛顿-柯特斯系数公式c(n,k)的详细推导。
当n>8时牛顿-柯特斯公式并没有实际意义。实际上，我们通常只用到n=1,2,4的情形，相应的公式分别称为梯形公式，辛卜生公式，和柯特斯公式。对于现代的计算工具来说，有梯形公式和辛卜生公式也就够用了。
利用牛顿-柯特斯系数，我们可以方便地写出牛顿-柯特斯求积公式：
4.梯形公式在牛顿-柯特斯求积公式中，如果我们取n=1,那么k可以取0和1。由此所形成的求积公式就是梯形公式。
由于得到的结果是梯形面积公式，所以称梯形公式。
辛卜生公式在牛顿-柯特斯求积公式中，如果我们取n=1,那么k可以取0，1，2，由此所形成的求积公式就是辛卜生公式。
柯特斯公式作为课外作业，大家可以取n=4,相应地k可以取0，1,2，3和4，仿照上面的方式，可以得到：
从而可进一步写出相应的求积公式，这就是柯特斯公式。
在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基于牛顿-柯特斯公式产生的序列，而龙贝格序列则不是，所以柯特斯公式还是很有意义的。
6.3复化中点公式法一般说来，在一个较大的积分区间上利用较高阶的牛顿-柯特斯公式虽然可以得到较高的代数精度，但实际效果并不好，道理也不难理解。
我们可以把积分区间划分为若干等分，在每个子区间上利用较低阶的牛顿-柯特斯公式，并把这些积分值相加。按照这样的思路所得到的求积公式我们统称为复化型求积公式。
复化中点公式复化中点公式也许最不为人们所注意，以至在一般的教科书中还没有这个名称，我们在后面将会看到，对于求数值积分来说，它实际上是最有用的公式。
把积分区间[a,b]划分为n等分，记x0,x1,…,xn为等分点，记[xj-1,xj]为第j个子区间,zj为区间的中点，j=1,2,…,n，记h=(b-a)/n,记Mn为所有子区间上利用中点公式所求得的积分值的和，那么我们有
（1）
2.复化中点公式的等价形式我们也可以把积分区间[a,b]划分为2n等分，记y0,y1,…,y2n-1,y2n为等分点，那么我们也可以把下标为奇数的点y2k-1,看成是区间[y2(k-1),y2k]的中点，k=1,2,…,n。(如图)
y0 y1 y2 y3 y4 y2(n-1) y2n-1 y2n
x0 x1 x2 xn-1 xn
这样，我们也可以把复化中点公式写成
（2）
提示：这两种表示方法我们在后面都会用到，所以一定要记住。不管公式的形式如何，编写求Mn的程序还是相同的。
3.复化中点公式(c语言代码)
在实际应用中，n通常取2的某个整数幂2k，我们可以把计算结果存放在数组M[N]中，程序段命名为GetM(int k),C语言源程序文件可以按如下的方式编写。
#include <math.h>
#define F(X) 4.0/(1.0+X*X)
static float a=0.0,b=1.0;
double GetMK(int k)
{ double x,y,step;
int n,j;
n=1;
for(j=0;j<k;j++) n*=2;
step = (b-a)/n;
x=a+step/2;
for(j=0;j<n;j++) { y+= F(x);x+=step;}
return (step*y);
}
6.4 变步长复化梯形公式法
1.复化梯形公式把积分区间[a,b]划分为n等分，记x0,x1,…,xn为等分点，记h=(b-a)/n,则有
记Ij,j=1,2,…,n为f(x)在第j个子区间[xj-1,xj]应用梯形公式所求得的积分值，则有
这就是我们的复化梯形公式。虽然我们也可以直接编写求Tn的计算机程序，但是没有这个必要。
变步长复化梯形公式假设对某个n,我们利用复化梯形公式，也就是上面的（3）式，得到了Tn，如果它不满足我们的精度要求，那么我们可以把每个子区间再对分一次，这相当于把积分区间划分为2n等分。
记y0,y1,y2,…,y2(n-1),y2n-1,y2n为等分点，记t=(b-a)/(2n),则有
再利用（3）式即得
亦即
（4）
这就是变步长复化梯形公式。我们可以利用它形成一个自动调整精度的算法。
变步长复化梯形公式方法为了方便地利用（4）式形成一个算法，我们总是取n等于2某个指数幂2k,并把Tn记为T[k],由n=2k可得2n=2k+1,从而有
T[k+1]=(T[K]+M[K])/2
算法说明如下：
T[0]=(b-a)*(f(a)+f(b)/2.0;
M[0]=GetMK(0)
K=1;
While(1)
{ M[k]=(M[k-1]+T[k-1])/2.0;
if(fabs(M[k]-M[k-1])<EPS)break;
GetMk(k);
K++;
}
return;
}
6.5 变步长复化辛卜生公式法
1.复化辛卜生公式假设x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-10,x2n-1,x2n把积分区间[a,b]划分为2n等分，我们也可以认为是其中的x0,x2,…,x2(k-1),x2k,…,x2(n-1),x2n把区间划[a,b]划分为n等份，并且x2k-1就是第k个子区间[x2(k-1),x2k]的中点。
记Ik,k=1,2,…,n为f(x)在第k个子区间[x2(k-1),x2k] 应用辛卜生公式所求得的积分值，则有
这就是我们的复化辛卜生公式,虽然我们也可以直接编写求S2n的计算机程序，但是没有这个必要。
2.关于复化辛卜生公式的结论假设x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-10,x2n-1,x2n把积分区间[a,b]划分为2n等分，那么我们可以得到下面3个不同的积分值：
利用x2k,k=0,1,…,n这n+1个点处的函数值和复化梯形公式计算出Tn ;
利用xj,j=0,1,…,2n这2n+1个点处的函数值和复化梯形公式计算出T2n ;
利用复化辛卜生公式计算出S2n ;
重要结论：对于上面的约定，我们有
S2n=(4T2n-Tn)/3
推导过程如下
4.变步长复化辛卜生公式加速算法我们仍然取n为2的整数幂的形式，即n=2k,并记Sn为S[k],那么我们有
T[k+1]=(M[k]+T[k])/2
S[k+1]=(4T[K+1]-T[K])
我们可以简单地对上一节中的变步长复化梯形公式法加以改进，几乎不增加每一步的计算量，即可得到一个加速算法。算法说明如下：
§5 龙贝格求积算法
1.统一的记号与公式对于求f(x)在区间[a,b]上的定积分值，在龙贝格求积算法中，我们采用统一的记号
它们之间的地推关系为
2.计算顺序对于上面给出的外推序列，我们可按下图所示的顺序进行计算，
如果有某个即可停止进一步的计算。
3.龙贝格算法的实现方法在龙贝格算法中，我们实际的做法是把梯形序列外推3次，即m值取0，1，2，3，并把相应的序列称为梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，并相应地记为T[k],S[k],C[k],R[k]。为了使我们的算法更为简明，我们再增加一个中点序列，记为M[k],并且用一个专门的函数GetM(int k)来计算（返回）M[k].
事实上，我们前面的变步长复化型辛卜生公式加速算法已经得到了中点序列，梯形序列和辛卜生序列。所以稍加改动即可得到龙贝格算法。算法说明如下
(2)
6.6 数值微分
1.利用差商替代微商求数值微分最简单的方法，也是最基本的方法，就是利用差商替代微商，余下的问题是如何估算误差。
假设f(x)是[x0-h,x0+h]上连续可微的实函数，当h的值充分小时，我们用f(x)在x0,x0+h处的一阶差商作为f’(x0)的近似值，即
f ’(x0)≈f[x0,x0+h]=[f(x0+h)-f(x0)]/h （1）
把f(x0+h)在x0处一阶泰勒展开后即可发现
|f ’(x0)-f[x0,x0+h]|=O(h)
此时我们称公式（1）具有一阶精度。
2．中点公式假如f(x)在x 0处适当光滑，也就是说我们所要用到的各阶导数都存在，那么我们可以把f(x0+h)和f(x0-h)在x0处泰勒展开得：
上面的（3）称为中点公式，也可以认为是利用中心差商替代微商。显然它具有二阶精度。
3.变步长方法我们也可以按变步长求数值积分的办法来求数值微分。假如对于给定的h>0,我们利用(3)式求f ’(x0)的近似值，并记为D[0],如果D[0]不满足精度要求，我们可以把h折半，重新进行计算。一般地，记
（4）
那么当K→∞时，我们有，D[k] →f ’(x0),如果我们利用某个|D[k+1]-D[k]|<EPS
作为停机规则，我们也可以得到一个自动调整精度的算法。
4．变步长加速方法对于求数值微分来说，我们可以简单地推出一个加速收敛的迭代格式，不过基本原理与上一节中的龙贝格求计算法相同，我们留到下一节中讨论。
在上面的（2）中分别取h为h/2k和h/2k+2，我们有
从而得到一个具有4阶精度的公式。
5.注意事项对于求数值微分来说，尽管从理论上讲步长愈小则结果愈好，但由于计算机的字长为固定值，所以当步长太小时会有下面一些潜在问题：
计算x0+h和x0-h时，可能出现大数吃小数的现象；
计算f(x0+h)-f(x0-h)时,可能出现绝对值相近符号相同的两个数相减，从而吃掉部分有效数字当分母太小时，商的绝对误差会大大增大。
总的说来，在当前的计算机条件下，这些问题的负面影响远没有想象的那么严重。
6.7里查逊外推法问题的提出假如我们要计算的某个真值A0可以表为与步长h有关的近似值形式F(t,,h),其中t是固定常数，且-H≤h≤H,h≠0或0<h≤H。
假设F(t,h)在h=0处可展为泰勒级数
F(t,h)=A0+A1h+A2h2+…+Aphp+… （1）
如果A1≠0,则称F(t,h)具有1阶精度；
如果A1=A2=…=Ap-1=0,Ap≠0,则称F(t,h)具有p阶精度；
结论：显然，一个近似计算公式F(t,h)的精度愈高愈好。
问题：我们能否利用利用低阶精度的公式构造高阶精度的公式？
2．中点公式法假如F(t,h)具有一阶精度或具有奇数p阶精度，而且允许h取负值，那么我们有
F(t,h)=A0+A1h+A2h2+…+Aphp+AP+1hp+1+… (2)
F(t,-h)=A0-A1h+A2h2+…-Aphp+AP+1hp+1+… (3)
从而进一步有
[F(t,h)-F(t,-h)]/2.0=A0+A2h2+…+AP+1hp+1+…
从而得到具有2阶或p+1阶的近似公式。
结论：如果一个近似公式具有奇数p阶精度，而步长又可以取负值，那么直接应用中点公式即可的到具有高阶精度的近似公式。
3．折半计算方法假如F(t,h)是一个具有p阶精度的近似公式，即
F(t,h)=A0+Aphp+O(hp+1) (4)
且不允许h取负值，那么我们可以用h/2替代h,从而有
F(t,h/2)=A0+Aphp/2p+O(hp+1) (5)
[(5)·2p-（4）]/(2p-1)得
[2p F(t,h/2)- F(t,h)]/ (2p-1)=A0+O(hp+1) (6)
结论：只要所给的近似公式适当光滑，我们总可以通过不同步长的线性组合来产生具有较高精度的近似公式。
术语：利用低阶公式的不同步长的线性组合来产生具有较高精度的近似公式称为外推。
4．里查逊外推法：
假如最初给出的近似公式F(t,h)只是1阶公式,-H≤h≤H,那么简单地利用中点公式即得：
[F(t,h)-F(t,-h)]/2.0=A0+A2h2+A4h4+…
记R0[k]= [F(t,h/2k)-F(t,-h/2k)]/2.0,k=0,1,…
R1[k]=[4·R0[k]-R0[k-1]]/3.0 k=1,2,…
R2[k]=(16·R1[k]-R1[k-1])/15.0 k=2,3,…
Rm[k]=(4mRm-1[k]- Rm-1[k-1])/（4m-1） k=m,m+1,…
则R0，R1，R2，…,Rm等序列的精度分别为2，4，6，….
5.求积公式的误差我们不必按照误差分析公式取确定各种复化型求积公式的误差，只需比较龙贝格算法和理查逊外推法的迭代格式即可记住：
复化梯形公式的误差为O(h2);
复化辛卜生公式的误差为O(h4);
复化柯特斯公式的误差为O(h6);
而龙贝格序列的误差为O(h8)。
6.1 求积公式由定积分的定义可知，连续函数f(x) 在区间[a,b]上的定积分近似值可以表示为[a,b]内的一些点x0,x1,…,xn处的函数值f(x0),f(x1),…,f(xn)的加权和或线性组合，即
（1）
其中w0,w1,…,wn仅与x0,x1,…,xn有关而与被积函数f(x)无关。我们把这样的公式称为求积公式，也称为机械求积公式。
1．术语和记号为了计算f(x) 在区间[a,b]上的定积分近似值，我们通常的做法是，把积分区间[a,b]划分为n等分，记h=(b-a)/n,x0=a,xk=a+kh,k=0,1,2,…,n,称x0,x1,…,xn为[a,b]的一个等份分划。
假如x0,x1,…,xn为[a,b]的一个等份分划那么求积公式(1)中的w0,w1,…,wn的选取仅仅只与n有关，从而可以简化对求积公式的研究。
2.求积公式的性质微积分学中我们曾研究过，定积分保持函数的线性关系不变，它的含义是，若f(x),g(x)都是[a,b]上的可积函数，则对任意实数u,v,我们有u·f(x)+v·g(x)也是[a,b]上的可积函数，而且
不难验证，求积公式也保持函数的线性关系不变，即
3.几种常见的求积公式在后面的讨论中，我们将经常用到下面一些非常简单的求积公式，他们是中点公式、梯形公式和辛卜生公式。
4．截断误差在求积公式中，我们使用的是近似等号，这是因为，对于一般的被积函数来说，利用这些公式计算所得的结果除了舍入误差外，还有截断误差。
有时为了进行误差分析，我们可以把求积公式式写成
（1’）
其中R[f]表示的就是截断误差。
5.代数精度的概念定义：一个求积公式
如果对所有的次数不超过m的多项式严格相等，而对某些m+1次多项式不相等，则称该公式具有代数精度m，或该公式的代数精度为m。
利用求积公式的线性性，我们不难证明下面的结论。
定理：如果求积公式对1，x,…,xm严格相等，而对xm+1不相等，则该公式的代数精度为m。
作为课外练习，鼓励大家给出完整证明。
6.基本结论我们不难利用上面的定理所给出的方法证明辛卜生公式的代数精度是3，而中点公式和梯形公式的代数精度是1。
现在我们可以对这三个公式作一个简单的评价：
中点公式和梯形公式的代数精度虽然都是1，但中点公式只计算一个点的函数值，而梯形公式要计算两个点处的函数值，所以中点公式优于梯形公式。
与梯形公式相比，辛卜生公式只多计算一个点的函数值，但代数精度却增加到3，显然辛卜生公式更为优越。
6．2 牛顿-柯特斯求积公式利用插值多项式近似替代被积函数设f(x)为被积函数，[a,b]为积分区间，x0,x1,…,xn为[a,b]内的n+1个互异的点，记Ln(x)为相应的拉格朗日插值多项式，那么我们有
2.利用插值多项式导出求积公式
在上面给出的公式中，由于诸lk(x)都是多项式函数，所以诸wk都可以精确地计算出来。从而我们可以得到一般性的求积公式。
3.讨论：
回顾上一章关于多项式插值的结论，由于任意次数不超过n的多项式与它的任意n+1个基点的插值多项式恒等，再由求积公式的代数精度的定义，我们立即得到：由n个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少式n。
为了同时说明求即公式得代数精度，我们记
4.牛顿-柯特斯求积公式牛顿-柯特斯求积公式就是利用等距基点的拉格朗日插值多项式导出的求积公式。
将积分区间[a,b]划分为n等分，记h=(b-a)/n，取
x0=a,xk=a+kh.k=0,1,…,n
我们可以得到
在牛顿-柯特斯公式中，我们称c(n,k)为牛顿-柯特斯系数，一般可通过查表得到。崔国华教材p60列出了直到n=8的所有牛顿-柯特斯系数，应该说，实用意义不大。崔国华教材p59-60给出了牛顿-柯特斯系数公式c(n,k)的详细推导。
当n>8时牛顿-柯特斯公式并没有实际意义。实际上，我们通常只用到n=1,2,4的情形，相应的公式分别称为梯形公式，辛卜生公式，和柯特斯公式。对于现代的计算工具来说，有梯形公式和辛卜生公式也就够用了。
利用牛顿-柯特斯系数，我们可以方便地写出牛顿-柯特斯求积公式：
4.梯形公式在牛顿-柯特斯求积公式中，如果我们取n=1,那么k可以取0和1。由此所形成的求积公式就是梯形公式。
由于得到的结果是梯形面积公式，所以称梯形公式。
辛卜生公式在牛顿-柯特斯求积公式中，如果我们取n=1,那么k可以取0，1，2，由此所形成的求积公式就是辛卜生公式。
柯特斯公式作为课外作业，大家可以取n=4,相应地k可以取0，1,2，3和4，仿照上面的方式，可以得到：
从而可进一步写出相应的求积公式，这就是柯特斯公式。
在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基于牛顿-柯特斯公式产生的序列，而龙贝格序列则不是，所以柯特斯公式还是很有意义的。
6.3复化中点公式法一般说来，在一个较大的积分区间上利用较高阶的牛顿-柯特斯公式虽然可以得到较高的代数精度，但实际效果并不好，道理也不难理解。
我们可以把积分区间划分为若干等分，在每个子区间上利用较低阶的牛顿-柯特斯公式，并把这些积分值相加。按照这样的思路所得到的求积公式我们统称为复化型求积公式。
复化中点公式复化中点公式也许最不为人们所注意，以至在一般的教科书中还没有这个名称，我们在后面将会看到，对于求数值积分来说，它实际上是最有用的公式。
把积分区间[a,b]划分为n等分，记x0,x1,…,xn为等分点，记[xj-1,xj]为第j个子区间,zj为区间的中点，j=1,2,…,n，记h=(b-a)/n,记Mn为所有子区间上利用中点公式所求得的积分值的和，那么我们有
（1）
2.复化中点公式的等价形式我们也可以把积分区间[a,b]划分为2n等分，记y0,y1,…,y2n-1,y2n为等分点，那么我们也可以把下标为奇数的点y2k-1,看成是区间[y2(k-1),y2k]的中点，k=1,2,…,n。(如图)
y0 y1 y2 y3 y4 y2(n-1) y2n-1 y2n
x0 x1 x2 xn-1 xn
这样，我们也可以把复化中点公式写成
（2）
提示：这两种表示方法我们在后面都会用到，所以一定要记住。不管公式的形式如何，编写求Mn的程序还是相同的。
3.复化中点公式(c语言代码)
在实际应用中，n通常取2的某个整数幂2k，我们可以把计算结果存放在数组M[N]中，程序段命名为GetM(int k),C语言源程序文件可以按如下的方式编写。
#include <math.h>
#define F(X) 4.0/(1.0+X*X)
static float a=0.0,b=1.0;
double GetMK(int k)
{ double x,y,step;
int n,j;
n=1;
for(j=0;j<k;j++) n*=2;
step = (b-a)/n;
x=a+step/2;
for(j=0;j<n;j++) { y+= F(x);x+=step;}
return (step*y);
}
6.4 变步长复化梯形公式法
1.复化梯形公式把积分区间[a,b]划分为n等分，记x0,x1,…,xn为等分点，记h=(b-a)/n,则有
记Ij,j=1,2,…,n为f(x)在第j个子区间[xj-1,xj]应用梯形公式所求得的积分值，则有
这就是我们的复化梯形公式。虽然我们也可以直接编写求Tn的计算机程序，但是没有这个必要。
变步长复化梯形公式假设对某个n,我们利用复化梯形公式，也就是上面的（3）式，得到了Tn，如果它不满足我们的精度要求，那么我们可以把每个子区间再对分一次，这相当于把积分区间划分为2n等分。
记y0,y1,y2,…,y2(n-1),y2n-1,y2n为等分点，记t=(b-a)/(2n),则有
再利用（3）式即得
亦即
（4）
这就是变步长复化梯形公式。我们可以利用它形成一个自动调整精度的算法。
变步长复化梯形公式方法为了方便地利用（4）式形成一个算法，我们总是取n等于2某个指数幂2k,并把Tn记为T[k],由n=2k可得2n=2k+1,从而有
T[k+1]=(T[K]+M[K])/2
算法说明如下：
T[0]=(b-a)*(f(a)+f(b)/2.0;
M[0]=GetMK(0)
K=1;
While(1)
{ M[k]=(M[k-1]+T[k-1])/2.0;
if(fabs(M[k]-M[k-1])<EPS)break;
GetMk(k);
K++;
}
return;
}
6.5 变步长复化辛卜生公式法
1.复化辛卜生公式假设x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-10,x2n-1,x2n把积分区间[a,b]划分为2n等分，我们也可以认为是其中的x0,x2,…,x2(k-1),x2k,…,x2(n-1),x2n把区间划[a,b]划分为n等份，并且x2k-1就是第k个子区间[x2(k-1),x2k]的中点。
记Ik,k=1,2,…,n为f(x)在第k个子区间[x2(k-1),x2k] 应用辛卜生公式所求得的积分值，则有
这就是我们的复化辛卜生公式,虽然我们也可以直接编写求S2n的计算机程序，但是没有这个必要。
2.关于复化辛卜生公式的结论假设x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-10,x2n-1,x2n把积分区间[a,b]划分为2n等分，那么我们可以得到下面3个不同的积分值：
利用x2k,k=0,1,…,n这n+1个点处的函数值和复化梯形公式计算出Tn ;
利用xj,j=0,1,…,2n这2n+1个点处的函数值和复化梯形公式计算出T2n ;
利用复化辛卜生公式计算出S2n ;
重要结论：对于上面的约定，我们有
S2n=(4T2n-Tn)/3
推导过程如下
4.变步长复化辛卜生公式加速算法我们仍然取n为2的整数幂的形式，即n=2k,并记Sn为S[k],那么我们有
T[k+1]=(M[k]+T[k])/2
S[k+1]=(4T[K+1]-T[K])
我们可以简单地对上一节中的变步长复化梯形公式法加以改进，几乎不增加每一步的计算量，即可得到一个加速算法。算法说明如下：
§5 龙贝格求积算法
1.统一的记号与公式对于求f(x)在区间[a,b]上的定积分值，在龙贝格求积算法中，我们采用统一的记号
它们之间的地推关系为
2.计算顺序对于上面给出的外推序列，我们可按下图所示的顺序进行计算，
如果有某个即可停止进一步的计算。
3.龙贝格算法的实现方法在龙贝格算法中，我们实际的做法是把梯形序列外推3次，即m值取0，1，2，3，并把相应的序列称为梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，并相应地记为T[k],S[k],C[k],R[k]。为了使我们的算法更为简明，我们再增加一个中点序列，记为M[k],并且用一个专门的函数GetM(int k)来计算（返回）M[k].
事实上，我们前面的变步长复化型辛卜生公式加速算法已经得到了中点序列，梯形序列和辛卜生序列。所以稍加改动即可得到龙贝格算法。算法说明如下
(2)
6.6 数值微分
1.利用差商替代微商求数值微分最简单的方法，也是最基本的方法，就是利用差商替代微商，余下的问题是如何估算误差。
假设f(x)是[x0-h,x0+h]上连续可微的实函数，当h的值充分小时，我们用f(x)在x0,x0+h处的一阶差商作为f’(x0)的近似值，即
f ’(x0)≈f[x0,x0+h]=[f(x0+h)-f(x0)]/h （1）
把f(x0+h)在x0处一阶泰勒展开后即可发现
|f ’(x0)-f[x0,x0+h]|=O(h)
此时我们称公式（1）具有一阶精度。
2．中点公式假如f(x)在x 0处适当光滑，也就是说我们所要用到的各阶导数都存在，那么我们可以把f(x0+h)和f(x0-h)在x0处泰勒展开得：
上面的（3）称为中点公式，也可以认为是利用中心差商替代微商。显然它具有二阶精度。
3.变步长方法我们也可以按变步长求数值积分的办法来求数值微分。假如对于给定的h>0,我们利用(3)式求f ’(x0)的近似值，并记为D[0],如果D[0]不满足精度要求，我们可以把h折半，重新进行计算。一般地，记
（4）
那么当K→∞时，我们有，D[k] →f ’(x0),如果我们利用某个|D[k+1]-D[k]|<EPS
作为停机规则，我们也可以得到一个自动调整精度的算法。
4．变步长加速方法对于求数值微分来说，我们可以简单地推出一个加速收敛的迭代格式，不过基本原理与上一节中的龙贝格求计算法相同，我们留到下一节中讨论。
在上面的（2）中分别取h为h/2k和h/2k+2，我们有
从而得到一个具有4阶精度的公式。
5.注意事项对于求数值微分来说，尽管从理论上讲步长愈小则结果愈好，但由于计算机的字长为固定值，所以当步长太小时会有下面一些潜在问题：
计算x0+h和x0-h时，可能出现大数吃小数的现象；
计算f(x0+h)-f(x0-h)时,可能出现绝对值相近符号相同的两个数相减，从而吃掉部分有效数字当分母太小时，商的绝对误差会大大增大。
总的说来，在当前的计算机条件下，这些问题的负面影响远没有想象的那么严重。
6.7里查逊外推法问题的提出假如我们要计算的某个真值A0可以表为与步长h有关的近似值形式F(t,,h),其中t是固定常数，且-H≤h≤H,h≠0或0<h≤H。
假设F(t,h)在h=0处可展为泰勒级数
F(t,h)=A0+A1h+A2h2+…+Aphp+… （1）
如果A1≠0,则称F(t,h)具有1阶精度；
如果A1=A2=…=Ap-1=0,Ap≠0,则称F(t,h)具有p阶精度；
结论：显然，一个近似计算公式F(t,h)的精度愈高愈好。
问题：我们能否利用利用低阶精度的公式构造高阶精度的公式？
2．中点公式法假如F(t,h)具有一阶精度或具有奇数p阶精度，而且允许h取负值，那么我们有
F(t,h)=A0+A1h+A2h2+…+Aphp+AP+1hp+1+… (2)
F(t,-h)=A0-A1h+A2h2+…-Aphp+AP+1hp+1+… (3)
从而进一步有
[F(t,h)-F(t,-h)]/2.0=A0+A2h2+…+AP+1hp+1+…
从而得到具有2阶或p+1阶的近似公式。
结论：如果一个近似公式具有奇数p阶精度，而步长又可以取负值，那么直接应用中点公式即可的到具有高阶精度的近似公式。
3．折半计算方法假如F(t,h)是一个具有p阶精度的近似公式，即
F(t,h)=A0+Aphp+O(hp+1) (4)
且不允许h取负值，那么我们可以用h/2替代h,从而有
F(t,h/2)=A0+Aphp/2p+O(hp+1) (5)
[(5)·2p-（4）]/(2p-1)得
[2p F(t,h/2)- F(t,h)]/ (2p-1)=A0+O(hp+1) (6)
结论：只要所给的近似公式适当光滑，我们总可以通过不同步长的线性组合来产生具有较高精度的近似公式。
术语：利用低阶公式的不同步长的线性组合来产生具有较高精度的近似公式称为外推。
4．里查逊外推法：
假如最初给出的近似公式F(t,h)只是1阶公式,-H≤h≤H,那么简单地利用中点公式即得：
[F(t,h)-F(t,-h)]/2.0=A0+A2h2+A4h4+…
记R0[k]= [F(t,h/2k)-F(t,-h/2k)]/2.0,k=0,1,…
R1[k]=[4·R0[k]-R0[k-1]]/3.0 k=1,2,…
R2[k]=(16·R1[k]-R1[k-1])/15.0 k=2,3,…
Rm[k]=(4mRm-1[k]- Rm-1[k-1])/（4m-1） k=m,m+1,…
则R0，R1，R2，…,Rm等序列的精度分别为2，4，6，….
5.求积公式的误差我们不必按照误差分析公式取确定各种复化型求积公式的误差，只需比较龙贝格算法和理查逊外推法的迭代格式即可记住：
复化梯形公式的误差为O(h2);
复化辛卜生公式的误差为O(h4);
复化柯特斯公式的误差为O(h6);
而龙贝格序列的误差为O(h8)。
