龙贝格积分实验报告

二、Romberg积分法

1.变步长Romberg积分法的原理

复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且f''(x)和f(4)(x)在区间[a,b]上的上界M的估计是较为困难的。在实际计算常采用变步长的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg积分法的思想。

在步长的逐步分半过程中，要解决两个问题： 1. 在计算出TN后，如何计算T2N,即导出T2N和TN之间的递推公式； 2. 在计算出TN后，如何估计其误差，即算法的终止的准则是什么。

首先推导梯形值的递推公式，在计算TN时，需要计算N1个点处的函数值在计算出TN后，在计算T2N时，需将每个子区间再做二等分，共新增N个节点。为了避免重复计算，计算T2N时，将已计算的N1个点的数值保留下来，只计算新增N个节点处的值。为此，把T2N表示成两部分之和，即

T2N2N11h2N[f(a)f(b)2f(akh2N)]2k1N1N1h2N[f(a)f(b)2f(akh2N)2f(a(2k1)h2N)] 2k1k1N1N1hN[f(a)f(b)2f(akh2N)]h2Nf(a(2k1)h2N)]22k1k1由此得到梯形值递推公式

T2NN1TNh2Nf(a(2k1)h2N)2k11

因此

h1[f(a)f(b)],2

11h2h1,T2T1h2f(ah2)22h1ba,T1. . . .

. .

由复化梯形公式的截断误差有

ITNIT2Nba2''hNf(1),a1b12 ba2''h2Nf(2),a2b12若f''(x)变化不大时，即f''(1)f''(2)，则有

I4T2NTN1T2N(T2NTN)4132

1式（2）表明，用T2N作为定积分I的近似值，其误差大致为(T2NTN)，

3因此其终止条件为

T2NTN

其中是预先给定的精度。

2.Romberg积分公式

将上述方法不断推广下去，可以得到一个求积分的序列，而且这个序列很快收敛到所求的定积分。记

(0)(1)TNTN，将区间N等分的梯形值。TNSN，将区间N等分的Simpson

(3)(2)RN，将区间N等分的Romberg。 TNCN，将区间N等分的Cotes。TN由其可构造一个序列{TN(k)}，次序列称为Romberg序列，并满足如下递推关系：

T(0)1ba1(0)baNba(0)[f(a)f(b)],T2NTNf(a(2k1)), 222Nk12Nk1)(k1)4kT2(NTNT,k1,2,k41kN

以上递推公式就是Romberg积分递推公式。

3.Romberg积分程序

1. 置N1，精度要求，h1ba； 2. 计算T1(0)ba[f(a)f(b)]； 2. . . .

. .

3. 置h2NhN1(0)baNba(0)，并计算T2NTNf(a(2k1))； 22Nk12N24. 置MN,N2N,K1;

k1)(k1)4kT2(MTM5. 计算T； k41kM6. 若 M1，则转（7）；否则置MM，kk1转（5）； 27. 若T1(k)T1(k1)，则停止计算（输出T1(k)），否则转（3）。

4.Romberg积分法的应用

function [T,n] = romb(f,a,b,eps) double R;

if nargin<4,eps=1e-8; end

h=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b)); n=1;J=0;err=1; while (err>eps)

J=J+1;h=h/2;S=0; for i=1:n

x=a+h*(2*i-1); S=S+feval(f,x); end

R(J+1,1)=R(J,1)/2+h*S; for k=1:J

R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1); end

err=abs(R(J+1,J+1)-R(J+1,J)); n=2*n; end R;

T=R(J+1,J+1) End

其中输入项：f为被积函数，ab为积分区间的端点值，ep为积分精度；输出项：T是逐次积分表值，n是迭代次数，R是最后积分值。 4.1程序调用

可以将被积分函数编成函数文件，也可以直接使用联函数来表示被积分函数，示例如下:

>>f=inline('1/(1+x.^2)','x');

. . . .

. .

>> [T,n,R]=romb(f,2,9,1e-9)

运行后得出其迭代次数，最终积分结果以及龙贝格积分矩阵如表2-1所示， 迭代次数N=，最终的积分值R=0.3530.

表2-1 龙贝格积分矩阵

0.7427 0.4833 0.3905 0.3628 0.3555 0.3536 0.3531 0.0000 0.3969 0.3596 0.3536 0.3530 0.3530 0.3530 0.0000 0.0000 0.3571 0.3532 0.3530 0.3530 0.3530 0.0000 0.0000 0.0000 0.3532 0.3530 0.3530 0.3530 0.0000 0.0000 0.0000 0.0000 0.3530 0.3530 0.3530 0.0000 0.0000 0.0000 0.0000 0.0000 0.3530 0.3530 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.3530 3.课本例题求解

21ln(1x)ln(1x)sin(x)(1)dx,(2)dx,(3)dx,(4)dx

21+x1xxx00001111 当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=0.6931

表2-2 式1对应的龙贝格积分矩阵

0.7500 0.7083 0.6970 0.6941 0.6934 0.6932 0.0000 0.6944 0.6933 0.6932 0.6931 0.6931 0.0000 0.0000 0.6932 0.6931 0.6931 0.6931 0.0000 0.0000 0.0000 0.6931 0.6931 0.6931 0.0000 0.0000 0.0000 0.0000 0.6931 0.6931 0.0000 0.0000 0.0000 0.0000 0.0000 0.6931 2 当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=0.2722.

表2-3 式2对应的龙贝格积分矩阵

0.1733 0.2488 0.2665 0.2708 0.2718 0.2721 0.0000 0.2740 0.2723 0.2722 0.2722 0.2722 0.0000 0.0000 0.2722 0.2722 0.2722 0.2722 0.0000 0.0000 0.0000 0.2722 0.2722 0.2722 0.0000 0.0000 0.0000 0.0000 0.2722 0.2722 0.0000 0.0000 0.0000 0.0000 0.0000 0.2722 3对于积分1ln(1x)dx，由于积分下限0为其奇点，理论上无法进行数值积分，x0本题中近似取下限为1*10-9来进行计算。当迭代精度ep=1e-9的条件下，迭代次

数N=16，迭代结果R=0.2722.

. . . .

. .

表2-4 式3对应的龙贝格积分矩阵

0.8466 0.8288 0.8241 0.8229 0.8226 0.0000 0.8228 0.8225 0.8225 0.8225 0.0000 0.0000 0.8225 0.8225 0.8225 0.0000 0.0000 0.0000 0.8225 0.8225 0.0000 0.0000 0.0000 0.0000 0.8225 4.对于积分sin(x)dx，同样积分下限0为积分函数的奇点，理论上无法进行数x02值积分运算，本题中仍取积分下限近似为1*10-9进行计算。当迭代精度ep=1e-9

的条件下，迭代次数N=16，迭代结果R=1.3708.

表2-5 式4对应的龙贝格积分矩阵

1.28 1.3498 1.3655 1.3695 1.3704 0.0000 1.3713 1.3708 1.3708 1.3708 0.0000 0.0000 1.3708 1.3708 1.3708 0.0000 0.0000 0.0000 1.3708 1.3708 0.0000 0.0000 0.0000 0.0000 1.3708 . . . .