完整版龙贝格积分试验报告

1、:、Romberg积分法1.变步长Romberg积分法的原理复化求积方法对丁提升精度是行之有效的方法,但复化公式的一个主要缺点 在丁要事先估计出部长.假设步长过大,那么精度难丁保证；假设步长过小,那么计算量 乂不会太大.而用复化公式的截断误差来估计步长,其结果是步长往往过小,而且f''x和fx在区间a,b上的上界M的估计是较为困难的.在实际计算中通常采用变步长的方法,即把步长逐次分半也就是把步长二等分,直到到达某种精度为止,这种方法就是 Romberg积分法的思想.在步长的逐步分半过程中,要解决两个问题：1. 在计算出Tn后,如何计算T2N ,即导出T2N和Tn之间的递推公式；

2、2. 在计算出Tn后,如何估计其误差,即算法的终止的准那么是什么.首先推导梯形值的递推公式,在计算 Tn时,需要计算N 1个点处的函数值在计算出Tn后,在计算T2n时,需将每个子区间再做二等分,共新增 N个节点.为了预防重复计算,计算T2N时,将已计算的N 1个点的数值保存下来,只计算新增N个节点处的值.为此,把T2N表示成两局部之和,即1,-h2Nf(a) f(b) 212h2Nf (a)1 hNNf(a)2 2T2Nf(b)f (b)N2 k N22N 12k 11f (a11f (af (a kh2N )%)kh2N)h2N由此得到梯形值递推公式T2N7Tn2h2N因此f(a (2 k

3、1膈)Nf(a (2 k 1)h2N) k 1Nf (a (2k 1)h2N)L L L L L L L L 1 k 1h1b a,T1.1 .十h2加,丁22#f(a) f(b), 2J ,、 T1 h2 f (a h2) 2由复化梯形公式的截断误差有b a 2 r" x,I Tn -hN f ( i), a 1 b 12 b a 21T2N. _ h2N f ( 2),a 2 b12假设f''(X)变化不大时,即f''( 1)f ''( 2),那么有I4T2N TnT2N 】(T2N Tn)L LLLLLLL24 13式(2)说明,

4、用T2N作为定积分I的近似值,其误差大致为-(T2N Tn), 3因此其终止条件为T2N Tn其中 是预先给定的精度.2.Romberg积分公式将上述方法不断推广下去,可以得到一个求积分的序列,而且这个序列很快 收敛到所求的定积分.记TN0) Tn ,将区间N等分的梯形值.Tf) Sn ,将区间N等分的Simpsont/2) Cn ,将区间N等分的Cote麟Tn Rn ,将区间N等分的Romberg.由其可构造一个序列TNk),次序列称为Romberg序列,并满足如下递推关系：铲¥f(a)晌福押)%a Nf(ak 1(2 kk (k 1) (k 1)T k 4 T 2 N T N k

5、1,2,Lk41以上递推公式就是Romberg积分递推公式.3.Romberg积分程序1. 置N 1,精度要求,h b a ；2. 计算 T1(0) Af(a) f(b);23.置h2NhN项计算试2tN0) Nb af (a2N k 14.N, N 2N, K 1;5.计算tM4kT2(M1)TMk1).416.7.假设 T1(k)M , ,k k 2,那么停止计算输出T1k,1 ,那么转7;否那么置M1 转(5);否那么转3.4.Romberg积分法的应用function T,n = romb(f,a,b,eps)double R;if nargin<4,eps=1e-8;endh=

6、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:nx=a+h*(2*i-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;for k=1:JR(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k)/(4Ak-1);enderr=abs(R(J+1,J+1)-R(J+1,J);n=2*n;endR;T=R(J+1,J+1)End其中输入项：f为被积函数,ab为积分区间的端点值,ep为积分精度；输出项：T 是逐次积

7、分表值,n是迭代次数,R是最后积分值.4.1程序调用可以将被积分函数编成函数文件,也可以直接使用内联函数来表示被积分函数,例如如下:>>f=inline('1/(1+x.A2)','x');>> T,n,R=romb(f,2,9,1e-9)运行后得出其迭代次数,最终积分结果以及龙贝格积分矩阵如表2-1所示,迭代次数N=64,最终的积分值R=0.3530.表2-1龙贝格积分矩阵0.74270.00000.00000.00000.00000.00000.00000.48330.39690.000010.00000.00000.00000.00

8、000.39050.35960.35710.00000.00000.00000.00000.36280.35360.353210.35320.00000.00000.00000.35550.35300.353010.35300.35300.00000.00000.35360.35300.353010.35300.35300.35300.00000.35310.35300.35300.35300.35300.35300.35303.课本例题求解1(1) : dx,(2)°1+x1ln(1 x)厂 dx,(3) x1)dx,(4)*dx1当迭代精度ep=1e-9的条件下,迭代次数N=32

9、,迭代结果R=0.6931表2-2式1对应的龙贝格积分矩阵0.75000.00000.00000.00000.00000.00000.70830.69440.00000.00000.00000.00000.69700.69330.69320.00000.00000.00000.69410.69320.69310.69310.00000.00000.69340.69310.69310.69310.69310.00000.69320.69310.69310.69310.69310.69312当迭代精度ep=1e-9的条件下,迭代次数N=32,迭代结果R=0.2722.表2-3式2对应的龙贝格积分矩

10、阵0.17330.00000.00000.00000.00000.00000.24880.27400.00000.00000.00000.00000.26650.27230.27220.00000.00000.00000.27080.27220.27220.27220.00000.00000.27180.27220.27220.27220.27220.00000.27210.27220.27220.27220.27220.27223对丁积分ln(1 x)dx,由丁积分下限.为其奇点,理论上无法进行数值积分,0 x此题中近似取下限为1*10-9来进行计算.当迭代精度ep=1e-9的条件下,迭代次 数N=16,迭代结果 R=0.2722.0.84660.00000.00000.00000.00000.82880.82280.00000.00000.00000.82410.82250.82250.00000.00000.82290.82250.82250.82250.00000.82260.82250.82250.82250.8225表2-4式3对应的龙贝格积分矩阵2表2-5式4对应的龙贝格积分矩阵1.28540.00000.00000.00000.00001.34981.37130.00000.00000.00001.36551.37081.37080.000