Romberg积分及初值问题

1、1 5.3 Romberg算算法法 综合前几节的内容,我们知道梯形公式,Simpson公式,Cotes公式的代数精度分别为1次,3次和5次复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为2阶、4阶和6阶无论从代数精度还是收敛速度,复合梯形公式都是较差的有没有办法改善梯形公式呢?2一、复合梯形公式的递推化等份分割为的积分区间将定积分nbadxxfIba,)(njjhaxk, 1 ,0,nabh各节点为 )()(2)(211njjnbfxfafnabT复合梯形(Trapz)公式为则不变等份，而分割为如果将,/)(2,nabhnba)()(2)(2)(41021112bfxfxfafn

2、abTnjjnjjn-(1)-(2)3)()(2)(2)(41021112bfxfxfafnabTnjjnjjnhjahxxjj)21(2121其中102111)(24)()(2)(4njjnjjxfnabbfxfafnab1021)(221njjnxfnabT10)21(221njnhjafnabT-(3)10)2)12(221njnnabjafnabT4abhn 时，1则由(1)(2)(3)式,有)()(21bfafabT)21(22112hafabTT)0(0T)1(0T)1(0kTTn记12kn若,2 , 1kjhaxj12kabhhxxjj212112)21(kabja12kabja

3、kabja2)12(5因此(1)(2)(3)式可化为如下递推公式)()(2bfafab)0(0T120001)2)12(2)1(21)(kjkkabjafabkTkT,2 , 1k(4)-上式称为递推的梯形公式6三种公式之间的关系244 1nnnTTS222441nnnSSC)1()(4141)1(11kTkTkTmmmmm)0(0T)1(0T)0(1T)2(0T)1(1T)0(2TRomberg算法求解步骤（3，-3）（-1,2）（0,1）（1，-1）（3，-4）7二、外推加速公式由复合梯形公式的余项公式)(3122nnnTTTInnTTI31342可得nnjjnTxfnabTI31) )(

4、221(341021由(3)式1021)(6)(431njjnxfnabT812kn设)1(31)(3400kTkT102111)(6)(4 )(2)()(231njjnjjxfnabxfbfafnabI )(4)(2)()(6102111njjnjjxfxfbfafnabnS复合Simpson公式nnTTI31342令引入),1(1kT)1(1kT)1(31)(3400kTkTnS12kS-(5)-(6)9)(15122nnnSSSI因此由复合Simpson公式的余项可得nnSSI15115162)1(1kT12kS即)1(151)(151611kTkTnS)(1kTnS2当然)1(2kT)

5、1(151)(151611kTkT令nC自己证明-(6)nC-(7)10)1(2kT12kCnC-(8)即)(2kTnC2当然同样由复合Cotes公式的余项)(63122nnnCCCInnCCI63163642)1(631)(636422kTkT得)1(3kT令)1(631)(636422kTkT-(9)11)1(1kT)1(31)(3400kTkT)1(151)(151611kTkT)1(2kT)1(3kT)1(631)(636422kTkT)()(2bfafab)0(0T120001)2)12(2)1(21)(kjkkabjafabkTkT,2 , 1k外推加速公式以上整个过程称为Romb

6、erg算法将上述结果综合后12)1()(4141)1(11kTkTkTmmmmm其中外推加速公式可简化为-(9)0(0T)1(0T)0(1T)2(0T)1(1T)0(2T)3(0T)2(1T)1(2T)0(3T,2 , 1mm可以推广到并且,2 , 1kRomberg算法求解步骤13romberg.m例:20sinxdx计算积分前侧矩形公式 z1 = 0.99212545660563 z11 = 0.99212545660563后侧矩形公式 z2 = 1.00783341987358 z22 = 1.00783341987358梯形公式 z3 = 0.99997943823961Simpson

7、公式 z4 = 1.000000002016138阶simpson公式 z5 =1.00000000000000自选步长梯形公式 z6 = 0.99999921563419自选步长Simpson公式 z7 =1.00000051668471Romberg公式 z8 = 0.99999999999802Mote-Carlo算法 z9 = 0.99821071589516 -0.00787454339437 -0.00787454339437 0.00783341987358 0.00783341987358 -0.00002056176039 0.00000000201613 -0.000000

8、00000000 -0.00000078436581 0.00000051668471 -0.00000000000198 -0.00178928410484Jifenbijiao.m积分法积分值绝对误差14如何构造Romberg算法3 龙贝格(Romberg)积分方法 我们已经知道,当被积函数f(x)在区间a,b上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分 可将分点(即基点)加密,也就是将区间a,b细分,然后利用复合梯形公式或复合抛物线公式求积。( )baf x dx 若用Tm表示把a,b作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h=h/2,则T

9、2m与Tm之间有如下关系: 2112(21) mmmkTThf akh(528) 其中 另外,若用Sm表示把a,b分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h就有20122221201222212024222(4224)3(244442)3(222)3MmmmmmmmmhSffffffhffffffhfffff 从Tm的定义可得到关系式 22441mmmTTS(529) 我们再举一个计算上半单位圆面积的例子(它的准确面积为/2)。现用内接正多边形的逼近方法来计算。 如图5.6,图(a) 、 图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是

10、用梯形方法计算其面积,图(b)是用三角形方法计算其面积。 图 5.6 设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为 ( )201235111133524(1)350112sincossin2sin242112()()23! 25! 222()()23!25!222()()23!25!2kkkkkkkkkkkkkkkknnTnnnT同理 如果组合一下,就会得到更精确的结果,即 (1)( )( )00154441225!(2 )kkkkTTT同理 (2)(1)(1)00151 4441225!(2)kkkkTTT 再以类似方法组合得 2(1)( )( )11164411()2(2

11、 )kkkkTTTO 这样继续下去,其值越来越接近上半单位圆面积/2。这种方法可以用到计算定积分 ( )baf x dx 为了推广公式(529)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。 龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。 首先将a,b依次作20,21,22,等分,记,0,1,2,2iibahi 按复合梯形公式(520)算得的值相应地记为T(k)0(k=0,1,2,);把按式(529)算得的S2m依次记为T(k)1(k=0,1,2,崐),而这每一个S2m又理解为由T2m与Tm的线

12、性组合得到的改进值,即(1)( )( )0014,0,1,2,41kkkTTTk 我们可按照类似的方法继续进行改进,也即由S2m与Sm的线性组合得到改进值,依次记为T(k)2(k=0,1,2,),即 2(1)( )( )10224,0,1,2,41kkkTTTk 这样就可构造出一个数表 (5-30) 其中除第0列(即最左一列)的T(k)0是按复合梯形公式计算外,其余各列都按下述规则(对m) (1)( )( )111,2,4,0,1,41mkkkmmmmmTTTk(531) 递推地计算出来。箭头表示计算流程。其计算步骤为: (1)将区间a,b等分为20,用梯形公式计算T(0)0,即(0)0 (

13、)( )2baTf af b (2)将区间a,b等分为21,用梯形公式算出T(1)0,即(1)0 ( )2 ()( )42baabTf aff b 再由T(0)0,T(1)0根据公式(531)算出T(0)1,即(0)(0)(0)101441TTT若 T(0)1-T(0)0, (为预给的精度) 则停止计算;否则继续往下计算; (3)依次分别算出T(2)0,T(1)1,T(0)2,这一行地往下推算,每一行算完,就得验证T(0)m(m=1,2,)是否满足预给的精度,即若(0)(0)1mmTT则停止计算;否则继续进行下一行。为了便于在计算机上实现,可运用下列公式编制程序:1(0)02( )(1)001

14、(1)( )( )11 ( )( )21(21)2221,2,4,0,1,41kkkkkiikkkiiiibaTf af bbabaTTf aiiTTTk 例 4 计算积分12041Idxx精确到10-4。 解(0)0(1)(0)00(1)(0)(0)00111(1) (0)(1)(42)322111(2)( )3.122243.13333341TffTTfTTT(2)(1)00(2)(1)(1)0012(1)(0)(0)11221113(3) ( )( )3.131176244443.1415694143.14211841TTffTTTTTT(3)(2)00(3)(2)(2)0012(2)(

15、1)(1)11223(1)(0)(0)2233111357(4) ( )( )( )( )2888883.13898844314159441443.14158641TTffffTTTTTTTTT(4)(3)00(4)(3)(3)0012(3)(2)(2)11223(2)(1)(1)223311135(5) ()()()21616161679111315()()()()()16161616163.14094243.1415934143.1415934143.14159341TTffffffffTTTTTTTTT4(1)(0)(0)3344(0)34(0)12011020

16、43.141593410.0000070.0000543.1415931443.14159261TTTTTdxxIdxarctgxx于是 由于 实际上 简单的数值方法与基本概念简单的数值方法与基本概念 1. 简单简单欧拉法（欧拉法（Euler） 2后退的欧拉法后退的欧拉法 3梯形法梯形法 4改进改进EulerEuler法法 2、初值问题的数值解法、初值问题的数值解法单步法单步法1. 简单的欧拉简单的欧拉(Euler)方法方法考虑模型：00( , ) (1.1)() (1.2)yf x yaxby xy 在精度要求不高时通过欧拉方法的讨论欧拉方法欧拉方法2. 欧拉方法的导出欧拉方法的导出把区间a

17、,b分为n个小区间步长为要计算出解函数 y(x) 在一系列节点()iiyy x)/()( iiixaihhhban,一般取即等距节点处的近似值01 naxxxb1(-)iiihxxN N等分等分00( ,)(1.1)() (1.2)yfx yaxby xy 对微分方程(1.1)两端从1nnxx到进行积分11( ,( )nnnnxxxxy dxfx y xdx11()()( , ( )nnxnnxy xy xf x y x dx右端积分用右端积分用左矩形数值左矩形数值求积公式求积公式2()( )()()()2bagg x dxba g aba( )( , ( )g xf x y x令1)1(,

18、()(, () nnnnxxnnf xy xnnyyf xy xh得x0 x11()()()(,)nnnnnny xy xhy xyh f xy)1,., 0(),(1 niyxfhyyiiii亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 每步计算只用到1ny或用向前差商近似导数1()()()nnny xy xyxh100021111(,)(,) (,)nnnnyyhf xyyyhf xyyyhf xy依上述公式逐次计算可得：ny故也称Euler为单步法。公式右端只含有已知项所以又称为显格式的单步法。 例例1 用欧拉公式求解初值问题用欧拉公

19、式求解初值问题)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步长取步长h=0.1，欧拉公式的具体形式为，欧拉公式的具体形式为)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy依次计算下去，依次计算下去，部分计算结果部分计算结果见下表见下表. xy21 与准确解与准确解 相比，可看出欧拉公式的计算结相比，可看出欧拉公式的计算结果精度很差果精度很差.

20、xn 欧拉公式数值解欧拉公式数值解yn准确解准确解y(xn) 误差误差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 欧拉公式具有明显的几何意义欧拉公式具有明显的几何意义, , 就是就是用折线近似用折线近似代替方程的解曲线代替方程的解曲线，因而常称公式，因而常称公式(2.1)为为欧拉折线法欧拉折线法. .( )yy xxynx1nxn

21、p1np1np x 还可以通过几何直观来考察欧拉方法的精度还可以通过几何直观来考察欧拉方法的精度. .假假设设yn=y(xn), ,即顶点即顶点Pn落在积分曲线落在积分曲线y=y(x)上，那么，上，那么，按欧拉方法做出的折线按欧拉方法做出的折线PnPn+1便是便是y=y(x)过点过点Pn的切线的切线. .从图形上看从图形上看, ,这这样定出的顶点样定出的顶点Pn+1显著显著地偏离了原来的积分曲地偏离了原来的积分曲线，可见欧拉方法是线，可见欧拉方法是相相当粗糙当粗糙的的. .定义定义若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。Ri 的的主项主项/* leading term

22、*/ 欧拉法的局部截断误差：232()()hiyxO h欧拉法具有 1 阶精度。),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR 定义定义在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。5. 欧拉公式的改进欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy

23、)1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxyih 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 设用欧拉公式设用欧拉公式),() 0(1nnnny

24、xhfyy 给出迭代初值给出迭代初值 ，用它代入，用它代入(2.5)式的式的右端，使之转右端，使之转化为显式，直接计算得化为显式，直接计算得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反复进行，得如此反复进行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn6.梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(21

25、11 niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyx

26、fhyyiiii假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量