第16 页 共 52 页
此式说明，积分的精确值近似地等于先后两个梯形积分的线性组合。这个新
1的积分近似值等于梯形积分T2n加上一个修正值(T2n?Tn)（误差补偿值）。
3事实上将上式右端中的Tn及T2n的表达式代入后，可以直接验证出这正是我们所熟知的辛浦生积分值。即：
41 Sn?T2n?Tn
33显然，这个新值有较快的收敛速度。这也就是说，利用梯形积分的先后两次积分近似值，进行恰当的线性组合，可以得到更为精确的辛浦生积分值。因此，这种方法有叫做线性加速法。
类似地，我们研究辛浦生公式的加速问题。
由于辛浦生公式的截断误差与小区间长度的四次方成正比，即有：
I*?S2n11 * ?4?I?Sn216将上式移项后解出I*得： I*?161S2n?Sn 1515不难验证，右端正是柯特斯积分Cn，即： Cn?1611S2n?Sn?S2?(Sn2?Sn) n151515此式说明，利用辛浦生积分的先后两次积分的近似值，恰当地进行线性组合，可以得到更为精确的柯特斯积分。
同样，对柯特斯积分公式也可以进行线性加速。
由于，柯特斯公式的截断误差与小区间长度的六次方成正比，由此可得：
I*?C2n1?6 *I?Cn2解出I*得：
I*?641C2n?Cn 6363这就是我们要推证的Romberg积分公式，并记作Rn，则可得：
第17 页 共 52 页
Rn?6411C2n?Cn?C2n?(C2n?Cn) 636363我们自然会想到，对Romberg积分再进行线性加速，不就能得到更为精确的积分公式了吗？实际上，这种想法是不现实的，回答是否定的。原因是，Romberg积分的插值多项式已经是八次，它具有九次插值多项式的精度。依据插值多项式的理论，当插值多项式的次数超过九次时，即将出现不稳定情况。因而，用线性加速的方法也就到此为止了。
Romberg积分的优点不仅精度高，而且它可以从最简单的梯形积分开始，逐步加工成收敛速度快的Romberg积分。
2.2 加速收敛技巧——Richardson外推法
假设有一个量F*，用一个步长为h的函数F1(h)去逼近，F*与h无关，对于F1(h)的截断误差有估计式
R(F*)?F*?F1(h)?a1hp1?a2hp2???akhpk?? (2.2.1)
pk?pk?1???p2?p1?0,ai(i?1,2,?）都是与h无关的常数。也就是说，F1(h)
逼近F*的阶是hp1，现在提出的问题是能否通过(2.2.1)构造出一个新的序列，他逼近F*的阶要比hp1更高，如为hp2。
将(4.4.1)中的h用qh来代替，则有
F*?F1(qh)?a1(qh)p1?a2(qh)p2???ak(qh)pk??
现在用hp1乘(2.1.7)的两边后和上式相减，整理得
(1?qp1)F*?(F1(qh)?qp1F1(h))?a2(qp2?qp1)hp2???ak(qpk?qp1)hpk??
若1?qp1?0，用1?qp1除等式两边有
ak(qpk?qp1)pka2(qp2?qp1)p2(2)p2(2)pk?ah??a?? (2.2.2) h??h????2khp1p11?q1?q其中
第18 页 共 52 页
a(2)2ak(qpk?qp1)a2(qp2?qp1)(2)?,?ak?,? p1p11?q1?q都是与h无关的常数，令
F1(qh)?qp1F1(h) F2(h)?1?qp1那么F2(h)逼近F*的误差由(4.4.2)知道为hp2。一般地当我们选定q为满足
1?qpm?0的适当正数时，计算公式为
Fm?1(qh)?qpm?1Fm?1(h)Fm(h)?,m?1,2,? (2.2.3)
1?qpm?1用归纳法容易证明，由(4.4.3)得到的Fm(h)逼近F*的误差为
(m)pm(m)pm?1F*?Fm(h)?amh?am?? ?1h(m)(m),am这里am?1,?是与h无关的常数。这就是Richardson外推法。
在后续的介绍中，我们将在复化梯形求积公式的基础上，用上述的Richardson外推法，构造Romberg积分的算法。
3 算法设计及数值实验
下面我们将重点给出各类积分公式的算法及程序设计，用他们来解决一些问题。在以下的程序设计中，我们将用Mathematica语言实现。
本文所取的数值例子也是经过精心挑选的。首先，为了能够比较计算值与精确值，我们必须采用可以用别的办法算出积分的那类例子，譬如，可以通过不定积分，变换，无穷级数，特殊函数表等等积出的例子。虽说这可能使人感到过于简单，但实际上这些例子是经过挑选以使他们确能反映出各种特性的。
3.1 梯形公式的算法思想及程序设计
在这里我们主要介绍两种常用的梯形公式：定步长梯形公式以及变步长梯形
第19 页 共 52 页
公式。我们将通过选用不同的例题来体现它们各自的用处和特点。
3.1.1定步长复化梯形公式
定步长的积分公式通常用于对于积分区间的长度有一定要求，或者是用于一些表函数，它的基本思想就是：当步长给定之后，在每一个积分的小区间上面运用梯形公式求积分。
下面我们来看一个例题：
例题1：用定步长的复化梯形公式来计算以下积分：
14i?01?x2dx, ?0x2dx 要求将区间分为2（i=1,2,?10）份。 11源程序及运行后结果见附录2。我们将其结果及精确解列表比较如下：
表二：定步长复化梯形求积方法示例 n 2 4 8 16 32 64 128 256 512 1024 精确值
由上表我们明显地看到梯形求积公式的收敛速度非常的低。而且，如果对于
4f(x)??dx 01?x21f(x)??xdx 01123.1000 3.1311 3.1389 3.1409 3.1414 3.1415 3.1415 3.1415 3.1415 3.1415 0.6035 0.6432 0.6581 0.6635 0.6655 0.6662 0.6665 0.6666 0.6666 0.6666 3.14159265 0.66666667 第20 页 共 52 页
精确度有明确的要求的话，此方法用起来显然有很多不便之处。
3.1.2变步长的复化梯形公式
对于变步长的复化梯形公式我们用逐次分半法处理。 取n=1，计算T1
T1=(b-a)(
f(a)f(b)) ?22把区间[a,b]分割为二份，取n=2计算T2.
T2=
Tb?ab?af(a)f(b)(??f(x1))?1?f(x1) 22222把区间[a,b]四等分，取n=4计算T4. T4= =
b?af(a)f(b)(??f(x1)?f(x2)?f(x3)) 422T2b?a?(f(x1)?f(x3)) 24每次总是在前一次的基础上将区间对分，分点加密一倍。愿分点上的函数值无需重复计算，一般计算公式为：
Tnb?anb?af(a?(2i?1)) T2n=??22ni?12n由复化梯形求积公式的误差估计式可以看到，当区间n等分时截断误差为 R(f,Tn)=-(b?a)2hf''(?n), a??n?b 12这里?n的下标n时表示?与分法有关。当取2n等分时，截断误差表示为
R(f,T2n)=-(b?a)h2()f''(?2n), a??2n?b 122将上面两式相减得到：
T2n-Tn=-b?ah2()[4f''(?n)?f''(?2n)] 122当f''(?)在区间[a,b]上连续，并假定n充分大时f''(?n)?f''(?2n)，则有