龙贝格积分实验报告.doc

二、 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时，将已计算的N1个点的数值保留下来，只计算新增N 个节点处的值为此，把 T2N 表示成两部分之和，即由此得到梯形值递推公式因此由复化梯形公式的截断误差有若 f'' ( x) 变化不大时，即f '' (1)f '' (2)，则有式（ 2）表明，用T2 N作为定积分 I 的近似值，其误差大致为1 (T2NTN ) ，因此其终3止条件为其中 是预先给定的精度。

2.Romberg 积分公式将上述方法不断推广下去，可以得到一个求积分的序列，而且这个序列很快收敛到所求的定积分记TN(0)TN ， 将区间N 等分的梯形值TN(1)SN ，将区间N 等分的SimpsonTN(2)C N，将区间N 等分的Cotes TN(3)RN，将区间N 等分的Romberg由其可构造一个序列 { TN(k ) } ，次序列称为 Romberg 序列，并满足如下递推关系：以上递推公式就是 Romberg 积分递推公式3.Romberg 积分程序1.置 N1， 精度要求， h1ba ；2.计算 T1(0)ba[ f (a)f (b)] ；23.置 h2NhN， 并计算 T2(0)N1 TN(0)b aNf ( a(2 k 1) b a) ；222Nk12 N4.置M N,N 2N,K 1;k(k1)(k1)5. 计算 TMk4T2 MkTM；41M6. 若 M 1， 则转（ 7）；否则置 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;endh=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)=(4^k*R(J+1,k)-R(J,k))/(4^k-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是逐次积分表值， n是迭代次数， R是最后积分值。

4.1 程序调用可以将0.74270.00000.00000.00000.00000.00000.0000被积分0.48330.39690.00000.00000.00000.00000.0000函数编0.39050.35960.35710.00000.00000.00000.0000成函数0.36280.35360.35320.35320.00000.00000.0000文件，0.35550.35300.35300.35300.35300.00000.0000也可以0.35360.35300.35300.35300.35300.35300.0000直接使0.35310.35300.35300.35300.35300.35300.3530用内联函数来表示被积分函数，示例如下 :>>f=inline('1/(1+x.^2)','x');>> [T,n,R]=romb(f,2,9,1e-9)运行后得出其迭代次数，最终积分结果以及龙贝格积分矩阵如表 2-1 所示，迭代次数 N=64，最终的积分值 R=0.3530.表 2-1 龙贝格积分矩阵3.课本例题求解1 当迭代精度 ep=1e-9的条件下，迭代次数 N=32，迭代结果 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 对应的龙贝格积分矩阵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.272213 对于积分0ln(1 x)dx ，由于积分下限 0 为其奇点，理论上无法进行数值积分，本题中近x似取下限为 1*10 -9 来进行计算。

当迭代精度 ep=1e-9的条件下，迭代次数 N=16，迭代结果R=0.2722.表 2-4式 3 对应的龙贝格积分矩阵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.822524.对于积分0sin( x)dx ，同样积分下限 0 为积分函数的奇点，理论上无法进行数值积分运算，x本题中仍取积分下限近似为 1*10-9 进行计算当迭代精度 ep=1e-9的条件下，迭代次数 N=16，迭代结果 R=1.3708.表 2-5式 4 对应的龙贝格积分矩阵1.28540.00000.00000.00000.00001.34981.37130.00000.00000.00001.36551.37081.37080.00000.00001.36951.37081.37081.37080.00001.37041.37081.37081.37081.3708。