5.3龙贝格(Romberg)算法

1第6章数值积分与数值微分龙贝格(Romberg)算法2综合前几节的内容,我们知道梯形公式,Simpson公式,Cotes公式的代数精度分别为1次,3次和5次复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为2阶、4阶和6阶无论从代数精度还是收敛速度,复合梯形公式都是较差的有没有办法改善梯形公式呢?龙贝格(Romberg)算法3等份分割为的积分区间将定积分nbadxxfIba],[)(njjhaxk,,1,0,nabh各节点为])()(2)([211njjnbfxfafnabT复合梯形(Trapz)公式为则不变等份，而分割为如果将,/)(2],[nabhnba)]()(2)(2)([41021112bfxfxfafnabTnjjnjjn--------(1)--------(2)一、复合梯形公式的递推化4)]()(2)(2)([41021112bfxfxfafnabTnjjnjjnhjahxxjj)21(2121其中102111)(24)]()(2)([4njjnjjxfnabbfxfafnab1021)(221njjnxfnabT10))21((221njnhjafnabT--------(3)10)2)12((221njnnabjafnabT5abhn时，1则由(1)(2)(3)式,有)]()([21bfafabT)21(22112hafabTT)0(0T)1(0T)1(0kTTn记12kn若,2,1kjhaxj12kabhhxxjj212112)21(kabja12kabjakabja2)12(即，k=1，n=1，k=2，n=2，k=3，n=4，k=4，n=8，……n=20=1n=1n=2n=21=26因此(1)(2)(3)式可化为如下递推公式)]()([2bfafab)0(0T120001)2)12((2)1(21)(kjkkabjafabkTkT,2,1k(4)-------上式称为递推的梯形公式递推梯形公式加上一个控制精度,即可成为自动选取步长的复合梯形公式具体的方法请同学们完成思考kT212kT12kT7由复合梯形公式的余项公式)(3122nnnTTTInnTTI31342可得nnjjnTxfnabTI31))(221(341021由(3)式1021)(6)(431njjnxfnabT二、外推加速公式812kn设)1(31)(3400kTkT102111)(6)(4])(2)()((2[31njjnjjxfnabxfbfafnabI])(4)(2)()((6102111njjnjjxfxfbfafnabnS复合Simpson公式nnTTI31342令引入),1(1kT)1(1kT)1(31)(3400kTkTnS12kS--------(5)--------(6)9)(15122nnnSSSI因此由复合Simpson公式的余项可得nnSSI15115162)1(1kT12kS即)1(151)(151611kTkTnS)(1kTnS2当然)1(2kT)1(151)(151611kTkT令nC自己证明--------(6)nC--------(7)10)1(2kT12kCnC--------(8)即)(2kTnC2当然同样由复合Cotes公式的余项)(63122nnnCCCInnCCI63163642)1(631)(636422kTkT得)1(3kT令)1(631)(636422kTkT--------(9)11)1(1kT)1(31)(3400kTkT)1(151)(151611kTkT)1(2kT)1(3kT)1(631)(636422kTkT)]()([2bfafab)0(0T120001)2)12((2)1(21)(kjkkabjafabkTkT,2,1k外推加速公式以上整个过程称为龙贝格(Romberg)算法将上述结果综合后-----(10)12龙贝格算法可用下表表示：13在龙贝格算法（10）式中，T0(k)表示梯形公式T2k,k=0,1,2,…T1(k-1)表示辛甫生公式S2k－1,k=1,2,3,…T2(k-2)表示柯特斯公式C2k－2,k=2,3,4,…T3(k-3)表示龙贝格公式R2k-3,k=3,4,5,…14)]1()(4[141)1(11kTkTkTmmmmm其中外推加速公式可简化为--------(11))0(0T)1(0T)0(1T)2(0T)1(1T)0(2T)3(0T)2(1T)1(2T)0(3T,2,1mm可以推广到并且,2,1kRomberg算法的收敛阶高达m+1的两倍Romberg算法求解步骤Romberg算法的代数精度为m的两倍15例3用龙贝格方法计算积分解：按(10)式计算公式计算如下1617程序框图开始计算T1、T2、T4、T8计算S1、S2、S4计算C1、C2计算R1|R1-C2|ε否k=4结束计算321-2222SkkkkRCT、、、||2322kkCR否k=k+1输出积分值I是是18龙贝格算法源程序#includestdio.h#includemath.hmain(){doublef(doublex);doublet0fun(intk,doublet0[50]);//计算T0(k)的函数说明doublet1fun(intk,doublet0[50],doublet1[50]);//计算T1(k)的函数说明doublet2fun(intk,doublet1[50],doublet2[50]);//计算T2(k)的函数说明doublet3fun(intk,doublet2[50],doublet3[50]);//计算T3(k)的函数说明doublet0[50],t1[50],t2[50],t3[50],a,b;//a,b为积分区间，t0[50]数组装T0(k)intk;//t1[50]数组装T1(k),依此类推。a=0.0;b=1.0;t0[0]=(b-a)*(f(a)+f(b))/2.0;//计算T0(0),即计算T119printf(\n\n);for(k=1;k=3;k++)//计算T0(1)、T0(2)、T0(3)、T1(0)、T1(1)、T1(2){t0fun(k,t0);//即计算T2、T4、T8、S1、S2、S4t1fun(k,t0,t1);}t2fun(1,t1,t2);//计算T2(0)，即计算C1t2fun(2,t1,t2);//计算T2(1)，即计算C2t3fun(1,t2,t3);//计算T3(0)，即计算R1if(fabs(t3[0]-t2[1])0.000000001)//若|R1-C2|ε，则停止计算printf(“I=%12.9f\n”,t3[0]);//输出积分近似值elsefor(k=4;k10;k++)//若不符合精度要求，则继续计算{t0fun(k,t0);//计算T0(4)，即计算T16(当K=4时)20t1fun(k,t0,t1);//计算T1(3)，即计算S8(当K=4时)t2fun(k-1,t1,t2);//计算T2(2)，即计算C4(当K=4时)t3fun(k-2,t2,t3);//计算T3(1)，即计算R2(当K=4时)if(fabs(t3[k-3]-t2[k-2])0.000000001)//若，停止计算{printf(k=%3d,I=%12.9f\n,k,t3[k-3]);break;//输出积分近似值，退出循环}}}doublef(doublex)//计算的函数{doubley;y=sqrt(1.0+x*x);returny;}||2322kkCR21doublet0fun(intk,doublet0[50])//计算T0(k)的函数，即计算的函数{intj,n;doubles,nk1,nk,a,b;nk1=pow(2.0,k-1.0);//计算2k-1nk=pow(2.0,k);//计算2ka=0.0;b=1.0;n=(int)(nk1-1.0);//取(nk1-1)的整数部分s=0.0;for(j=0;j=n;j++)//计算{s=s+f(a+(2*j+1.0)*(b-a)/nk);}t0[k]=0.5*t0[k-1]+(b-a)*s/nk;//计算T0(k)returnt0[k];}kT2)2)12((1201kjkabjaf22doublet1fun(intk,doublet0[50],doublet1[50])//计算T1(k-1)的函数，{t1[k-1]=(4.0*t0[k]-t0[k-1])/3.0;//即计算的函数returnt1[k-1];}doublet2fun(intk,doublet1[50],doublet2[50])//计算T2(k-1)的函数，{t2[k-1]=(16.0*t1[k]-t1[k-1])/15.0;//即计算的函数returnt2[k-1];}doublet3fun(intk,doublet2[50],doublet3[50])//计算T3(k-1)的函数，{t3[k-1]=(64.0*t2[k]-t2[k-1])/63.0;//即计算的函数returnt3[k-1];}12kS12kC12kR