华中科技大学计算方法课件-第六章常微分方程初值问题数值解法.pdf

第六章常微分方程初值问题的数值解法考虑d维常微分方程初值问题瑞=/(力，沙(力)，力eQ,b丘=()，当右函数/(力g)e C(a,b x H)满足Lipschitz条件I I M z/1)-fy2)y(tn+i).(4)采用先前约定符号即得求解系统(1)的 线 性 方法yn+1=Vn+h 0 fn+(1-0)fn+i.(5)方法(5)在求当前值勤+1时，仅需要先前一步的信息机，力 2，因此我们称这种方法为单步法.特别，若e=1 则得显式E ul er法Un+l=Vn+h f n(6)若J =0,则得隐式E ul er法V n+l=Vn+h fn+1(7)若。=I，则得梯形法hV n+l=Vn+fn+l)(8)方法(6)是一个显式方法，即其每个计算步只需直接递推就可获得当前数值解枷+1，而方法(7)及(8)是隐式方法，其每个计算步均需解一个关于垢+1的隐式方程才可获得当前数值解.若记l,n=yn,丫2,n=yn+hef(tn,H,n)+(1 -：)/(5+.迄 n),则方法(5)也可等价地写成二级Runge-Kutta方法,n yi)k n=Vn+H,n)+。+%石 一),(9)Vn+1=yn+hOf(tn,H,n)+(1 -O)f(tn+h,K,n)j.例 6.1 试利用显式Euler法及梯形法构造出一个改进的计算格式，并分别用显式Euler法及所构造的新格式计算初值问题(取步长九二0.1)(%)=作沙+。+exp(川 t e 0,1 ,5(0)=1,3且比较两者的计算精度。解 利用显式Euler法及梯形法可构造出如下改进的Euler法Vn+l=Vn+Vn)+f(力 n+1,Vn+(1 D取步长h=0.1,并分别应用方法(6),(11)于初值问题(10)可得如下数值结果tn数值解反 差Euler 法改进的Euler法Euler 法改进的Euler法0.11.30001.33503.7257e-0022.2122e-0030.21.67011.75388.873 le-0024.9783e-0030.32.12432.27291.5694e-0018.3367e-0030.42.67932.91172.4471e-0011.2326e-0020.53.35443.69263.5521e-0011.6984e-0020.64.17264.64234.9199e-0012.2350e-0020.75.16075.79136.5907e-0012.8462e-0020.86.34987.17548.6097e-0013.5360e-0020.97.77648.83611.1028e+0004.3080e-0021.09.482910.82151.3903e+0005.1663e-002其中初值问题(10)的精确解=(力+1产exp，网格点tn=nh处的误差为ytn)-yno由上表可见，改进的Euler法的计算精度要优于显式Euler法。6.1.2数值积分法数值积分法的基本思想是先将问题转化为积分方程y(tm)y(tn)=/(力 贝 力)此 m n,(12)J h然后采用第五章介绍的数值积分公式将上式右端离散化，从而获得原初值问题的差分格式。取?7 1 =八+1,应用矩形公式于(12)中，则 得 单 支 方 法Vn+i=ynhfOtn-(l O)tn+xOyn-(l O)yn+i,0 E 0,1.(13)特别，当9=2时，称之为隐式中点公式。若记H,n=Vn+h(l 0)ftn+(1 0)h,Y 1,n,则方法(13)也可以写为单级Runge-Kutta方法JK,7 2 =物+九(1 0)ftn+(1 0)h,Y 1,n,yn+l=Vn+hftn+(1 H,n.在公式(12)中取m=n +k(k e N),并将Newton-Cotes公式代入其中，则可得一类线性k步公式Vn+k=Un+h):Bifn+i,八=0,1,2,，N k,(15)2=0其 中，优=易 吝 灯n (“加 特 别，当取k=2时，可j=U,/i得Milne公式Vn+2=Un +(/n+2+fn+1+fn)-(16)o利 用 公 式(15)计 算 当 前 步 枷+斓 寸，其 需 要 前 七 步 逼 近值%2,枷+1),，yn+k1的 有 关 消 息，故称之为k步公式。关于这类公式更广泛的形式将在下段继续讨论。6.1.3 Taylor展开法Taylor展开法的基本思想是首先构造一个关于真解及其有关信息的含参算子，将算子中诸项在某点处按Taylor展式展开并截去余项，然后令诸同类项系数为零，由此确定出原算子中的全部或部分参数，从而获得一个或一类关于数值解的差分格式。设期向为待定系数，定义算子kLy(t),h=Q渭(力 +ih)-hyt+ih).(17)分=0若V有0+2阶连续导数，则由Taylor展式有p川川=+&+1犷+1/+1)+。(从+2),(18)j=o其中f k k。0=T%=52(沁,；Pi)_ Bo,i=E (19)G=(沁。一鸡)，=2,3,，0+1.氏 于n+t,其中 Q/c#0,舄+而#0.(22)z=0 i=0若(20)成立，则称方法(22)为阶相容的，此时称式(21)的右端为该方法的局部截断误差.特别，若%#0,则称之为隐式的，否则称之为显式的。上述分析过程表明：定理6.1 线性k步法(22)为p阶 相 容 的 充 要 条 件 是(20)成立,且4+140.在(20)中取k=2,p=3,3=1 得(C VQ+QI+1=0,Q 1+2(仇+仇+防)=苏 3 1 +4)(B i+2 =0,、疝 3 +8)一东(为 +4/52)=0,解之得12 1因=-1 Q o,0o=一 乃(1+5。0),f3 i=-(1 QQ),优=行(5+。0).J L /O _ L /从而得二步方法h如+2(1+。0)%1+1+。0472=行(5+0)/计2+8(1a o)/n+i (l+5a()/J.(23)对于该方法，我们有。4-(1 +必)，。5=一熹(17+13a o).24 3o(J因 此 当 a。#1时，。4#0，此时方法(23)为三阶的。特别，若取仅)=5,则得二步显式方法lJn+2+4 枷+i 5yn=2/z(2/n_|-i +于n)；(24)若取Q o=0,则得三阶A d a m s-M on l ton 方法Vn+2 Vn+1=(5/n+2+8/n+l fn)；(25)J L乙若取小=1,则。4=0,但 踊=击#0，此时方法(23)即为四阶M i l n e方法Vn+2 Vn=(/n+2+fn+1+fn (26)类似地，若在(20)中取k=0=3,Q3 q2=1，因=小，贝 U得三步Adams方法hVn+3 Vn+2+-(5 12仇)九+3+(8+36(3。)fn+2_ L 乙-(1+36 仇仇+1+12%加，(27)且此时C4=仇 五，。5法(27)为三阶的；当仇=即为四阶Adams方法9.oOo豺-+C 4445，-时=124因此,但。5当优#0,*吉 时，方此时该公式yn+3 Vn+2+何9/九+3+19/n+2 fn+1+fn-(28)6.2 Runge-Kutta 方法一个k步方法当用于计算初值问题(1)时，除需问题本身的初值加外，还需要k-1个额外的计算启动值以,yi)一)yk i*而一般单步方法则可自起始完成计算 但是，如果单步方法的每步计算：Vn 一%+1仅由前一步的逼近值姗直接导出，则这样的单步方法往往计算精度不高.如：线性。-方法(5)与 单 支 方 法(13)的局部截断误差均为Rn=(。一 J的/陶)+。仇3),则当e*外寸,其方法的相容阶仅等于1；当。=外寸,其方法即分别为梯形法写隐式中点法,其相容阶达到这类方法的最高阶2.为提高单步法而精度，Runge与Kutta分别提出了在由外到yn+1的许算过程中增加若干中间逼近值的数值计算方案,这就形成了我们本节将要介绍的Runge-Kutta方法.值得注意的是：在上节，我们已指出 线 性 方 法 与 单 支。-方法可写成Runge-Kutta方法形式，但其属低阶方法,本节介绍的Runge-Kutta方法将涵盖高阶方法.6.2.1 方法的一般形式Runge-Kutta方法的一般形式为S yn h 于Un+Cjh,Y j,n),7=1sVn+1=Vn+h b j f (tn+Cjh,X,n),0,h 0为 步 长，力几=对于方法(29),若5 j 时均有%=0,则称之为显式方法，否则称之为是隐式方法.为简化Runge-Kutta方法的书写，我们经常采甬所谓Butcher表c I Ab(30)来表征方法(2 9),其中A=(%)G R-s xs,c=(cb。2,cs)r,b=(5i,如，b,s)，e 瞪.一个Runge-Kutta方法(29)称为是 P阶相容的(或简称为2阶的),若其局部离散误差dn+l-=y(tn+l)-yn+1=O(hp+1 九 一。(31)其中 歹九+1由 0(32)确定.Runge-Kutta方法的相容阶概念是构造具体方法的基本点.6.2.2显式方法显式Runge-Kutta方法与隐式Runge-Kutta方法相比较，由于前者在每步无须解s个隐式方程，因此其具有计算简单的特点.具体显式Runge-Kutta方法的推导可采用Taylor展并法获得.下面,我们以三级三阶显式Runge-Kutta法为例来说明其导出过程.设V为系统(1)的充分可微解，则由Taylor展式有3 1y(tn+九)=y(tn)+V+。仇 4)z=i z!(33)=y(tn)+hfn+-h2 Fn+-h3Fnf+Gn)+。(九4)2o其中=/(%ri,?/(%),fn=df(tn,y(tn)ay?f n Jn Jri)Gn=a2 K te y g)dt2,9；啊(tn,y(tn)荏 啊(5,y十dGtd4yJn 彻2以下恒设三级显式Runge-Kutta法满足ci=0,2-1G):。V,2=2,3.J=I(34)若方法自精确值式5)出发计算一步，则有fl,n=/(一,9),于2,n=f (%+C?h,y(tn)+6a21/1,n),于3,n=f 8 +c3h,y(tn)+6a31/1,72+32/2,n),(35)3yn+i=u(17 2)+h bjfj,n.、1=1将介加了3外在(G g陶)处展开，并注意到条件(34)得f f2,n=A +hc2Fn+却 26Gzi+。仇 3)h,n=/n+hcFn+%2(。2。32 工+1c3n)+。户).将上述所获力,九(i=1,2,3)的表达式代入(35)的最后一式得yn+i=y(tn)+Hbi+3 +b3)fn+hb2C 2+b3c3)Fn1 Q -,9 9 /(37)+5力2b3c2a32月I北+02。2 +3c3)n+。(万).乙由相容阶定义及式(33)与(3 7),若要方法有三阶，则必有(34)及F列条件成立：(38)瓦+向+后=1,b2 c2+6363=I，2363C+22C2b1-6b3c2as2-取(34),(38)的一组解13瓦=j,历=0,63=储C i =。31=,。2=。21=1235 C 3=Q 3 2一,可得三阶Heun方法OO9)oOO2-31-3CII)O1-32-33-4O1-4取(34),(38)的另一组解1 2 1bl=b3=62=Cl=0,C2=Q 21=C 3=-231=1,32=2,0 3 2得三阶Kutta方法0 0 0 0；.1 2 0(4。)-16 3 6仿上,我们还可获得如下四级四阶方法0 0 0 0 00000019190 0 01212000120120 012y/2-l2-200(41)1 0 0 1 01022+24 016131316162-V262+/616上二方法分别称为经典Runge-Kutta方法与Gill方法。当构造显式Runge-Kutta方法时，下述由Butcher给出的最大可达阶定理是值得注意的。定理6.2 一个6级显式Runge-Kutta方法的相容阶不可能超过s,且不存在五级五阶显式Runge-Kutta法。以下是Gill方法的计算实现。算 法6.1 Gill方法z=gill(a,b,y O,d,h)A=0 0 0 0;1/2 0 0 0;(sq rt(2)-1)/2 (2-sq rt(2)/2 0 0;.0-s q rt(2)/2 (2+s q r t(2)/2 0;c=0 1/2 1/2 1;B=l/6(2-sq rt(2)/6(2+sqrt(2)/6 1/6;t0=a;for n=a+h:h:b11=tO+c(1)*h;t2=t0+c(2)*h;t3=t0+c(3)*h;t4=t0+c(4)*h;Yl=y0;Y2=yO+h*kron(A(2,1),eye(d)*f(tl,Y1);Y3=y0+h*kron(A(3,1:2),eye(d)*f(11 ,Y1);f(t2,Y2);Y4=yO+h*kron(A(4,1:3),eye(d)*f(11 ,Y1);f(t2,Y2);f(t3,Y3);yl=yO+h*kron(B(l:4),eye(d)*f(tl,Y1);f(t2,Y2);f(t3,Y3);f(t4,Y4);tO=tO+h;yO=y1;end例 6.2 设置在某边防岛屿的雷达发现正东方向15海里处有一艘敌舰侵入我领海，其正以每小时4 0 海里的速度向正北方向行驶，驻岛部队立即派出一艘舰艇以每小时4 5 海里的速度追击。试建立我舰艇的航海轨迹模型，描绘出其轨迹图，并问经多少时间我舰艇能追上敌舰？解 以边防岛屿为原点。(0,0),敌舰初始位置记为4 点，线段0 4 所在直线为力轴，过原点。且垂直于z轴的直线为v 轴，建立直角座标系如图(见图6.1)。设 敌 我 二 舰 在 力 时 刻 的 位 置 分 别 为D,5。由于6(/4)点的运动方向总是指向敌舰位置。，因此，向 量 标的方向即为我舰运动轨迹的切向，由图6.1中直角三角形 6C O的边角关系可得我舰的航海轨迹模型=0.0111 =tO+c(1)*h;t2=t0+c(2)*h;t3=t0+c(3)*h;t4=t0+c(4)*h;Yl=y0;Y2=yO+h*kron(A(2,1),eye(d)*f(tl,Y1);Y3=yO+h*kron(A(3,1:2),eye(d)*f(tl,Y1);f(t2,Y2);Y4=yO+h*kron(A(4,1:3),eye(d)*f(tl,Y1);f(t2,Y2);f(t3,Y3);yl=yO+h*kron(B(1:4),eye(d)*f(tl,Y1);f(t2,Y2);f(t3,Y3);f(t4,Y4);dist=abs(y 1 (1)k);hold on;p lo t(y l(l),y l(2),5b-.5)hold o ff;x la b e l(x);y la b e l(y);T=t0+h;tO=tO+h;yO=y1;end图 6.2 Gill方法解初值问题(42)的解曲线图.这里，我们考虑到计算机的舍入误差，没有取程序终止准则为dist=0,而是将其设置为dist 0.01（“dist”为敌我二舰之距离）。采用命令giUex（0/0 0；2,0.001）即得我舰艇运动的轨迹图（见图6.2）,其追上敌舰的时间T 七0.9230小时。6.2.3隐式方法由上节可见，随着方法的级数与阶数的增加，利用Taylor展式来利用Tay lor展式来Runge-Kutta方法获得愈来愈加困难，不言而喻，直接利用Taylor展式来获取隐式方法将更加困难。为克服该困难，Butcher于1964年引入了下列简化条件：f b吐1=1 =i=ls Ic:=亍，z=1,2,s;Z =1,2,?7j=is。：94T =7%(1 4)，=1,2,，s；/=1,2,，您一 Vi=l这里，B(p),。)分别意味着系统的真解沙满足y(tn+九)=y(tn)+h bjytn+CjK)+7=1y(tn+Cih)=y g+物”(如 +Cjh)+。(a+1).j=i基于上述阶简化条件,Butcher给出了Runge-Kutta方法的相容阶判据。定理6.3 若Runge-Kutta方法(29)满足阶简化条件。)，。仁)，且9为满足石(夕)及不等式p W min +(+1,2 +2的最大正整数,则该方法为夕阶相容的.依据阶简化条件，人们将常用隐式Runge-Kutta方法分为六类.表6.1 六类隐式Runge-Kutta方法方 法阶条件相 容 阶GaussB(2s),C(s),D2sRadau IAB(2s-1),1),O(s)2s-1Radau IIAB(2s-1),C(s),D(s-1)2s-1Lobatto IIIAB(2s-2),C(s),D(s-2)2s-2Lobatto IIIBB(2s-2),C(s-2),D(s)2s 2Lobatto IIIC B(2s 2),C(s 1),D(s 1)2s 2利用阶简化条件及定理6.3,我们可构造隐式Runge-Kutta方法.如：对于二级Runge-Kutta方法,若有以3),。，。成立，即2 2 优d 1 =%/=1,2,3;djj=Q,i 1,2,3 17=1(43)LuE Q，jbi%(1 Cj),j=1,2,0为积分步长，k 1是方法的步数，诸a7为实常数，且供*。标=a+nh,Nh b-a,yn y g,夕/:t 如是依赖于问题的右函数/的映射，其中DH=(加 Vn,yn+L-，yn+k；h)0 h H,a tn 0是适当选取的常数。恒设1.当/=0时，(pf=0；2.当/满足(2)时，灯满足Lipschitz条件1夕/（力 九;Un+l、,Vn+ki 卜）一夕/（力 疝 Zn,砺+1,一，无）|k Lq|Hn+i n+i?。%S 九 夕）2=0（45）这里却和4通常依赖于(2)中的Lipschitz常数L理 论 分 析 表 明：只 要 九 mm 雪,&,则从任意初值？/o,勿,，%_1出发，方程(44)可唯一地阖定问题的一个数值解序列 勿 黑0.方法(44)是一类非常广泛的方法，它几乎包含当今各种常用的各种常微分方程数值方法。如取kD/=优/(5+/，第 1+)i=0则得线性多步法(22)；若取夕/=):bj f (tn+Cjh,Y j,n),j=i贝 IJ得s级Runge-Kutta方法(29).此外，(44)也包括6.4中将要介绍的预校算法及迄今为止尚待开发的各种方法。本节将涉及该类方法的相容性、稳定性及收敛性。6.3.1相容性6.1和6.2中局部截断误差及相容阶概念可统一定义如下：定义6.1设/满足(2),若问题的真解？/满足k(+?；)=力 九；U(九+1),。y(tn+k)；九)+Tn+k)i=0n=0,1)2)N (46)则称其余项乙+k为方法(44)在点力计卜处的局部 截 断 误 差。进一步，若夕为满足max Tn+k =O(hp+1),无-0(47)0nN k的最大正数，则称方法(44)是p阶相容的。特别若r max|乙+上|0,/i 一 0,(48)h QnN k则称方法(44)为相容的。定理6.4方法(44)相容的充要条件是p(l)=0,且下列等式关于八一致成立J i m 0/陶;g(a+i),，y(tn+k)；h)=p娟陶),(49)九一0k其中 P=Q (S3。2=0由定理6.4可直接推得推论6.1线性多步法(22)相容的充要条件是P =。=。，k其中=(e C).2=0S推论6.2 Runge-Kutta方法(29湘容的充要条件是%=1.j=i6.3.2稳定性定义6.2 方法(44)称为是零稳定的，若对任意满足Lipschitz条件的问题，存在正常数。，厩 0 便 当 尻 时,方 法(44)的任一解序列 枷与相应的扰动问题(设诸7 历任给扰动)eg/U n；Z g *,%n+k；九)+7 2 =0,1,-,N k,Vj+50 J=0,1,，上 一 1,的解序列 为 满足翳 瞌 一 斯|摩。0舞I闾卜上述稳定性概念充分刻画了当步长九充分小时计算过程中扰动对数值解的影响。显然，依据定义6.2来判断零稳定是困难的，下述结论可在一定程度上简化其判定。定理6.5 方法(44)为零稳定的充要条件是多项式p(S)的每个根的模不超过1,且模为1的根是单根。方法(44)满足上述充要条件也称其符合根条件。由于对于任何单步方法有p(。=-1,因此单步方法必满足根条件，即是零稳定的。进而可知，Runge-Kutta方法均是零稳定。定理6.6(Dahlquist)零稳定的线性k步法(22)的相容阶p满足(卜+2,若k是偶数；p OO为方法(44)的绝对稳定域。作为范例，考虑线性多步法(22)和Runge-Kutta方法(2 9),当其应用于(50)时，分别产生下列差分方程kk iVn+i=力 ：Bijn+i)(51)yn+1=1+hbTI 一 hA)reyn,(52)其中，B=h入/为s级单位阵，e=(1,1,.,I f G Rs.若设上述差分方程有形如珈=居+0)的解，将其代入则分别得其特征方程为P=W(53)=1+田 q 和 加 一 上.(54)由此可知，线性多步法(22)和Runge-Kutta方法(29)的绝对稳定域分别为LM=h E C 方程p(S)=的根模满足|冬 1,SRK=h e C|1+hbTI-hA)-1e 1.6.3.3相容性当以方法(44)求解满足Lipschitz条件的问题时，若产生唯一的逼近序列彳物,则 称 误 差n-=y(tn)yn,(n=0,1,N),为方法(44)的整体截断误差。定义6.5 方法(44)称为P 阶 收 敛 的，若以该方法求解任意满足Lipschitz条件的问题(1),当/充分连续可微且max则）%|=6（次）,h t 0时，有 金=。(游)，h T 0。特别，方法称为是收敛的，若V力Ca,b,有yn ，当八一 0,a+nh-t,仇一 纨)(0 S 分 S k 1).定理6.7 若方法(44)相容，则其零稳定性与收敛性等价。定理6.8 若方法(44%阶相容且零稳定，则该方法是夕阶收敛的。由推论6.1,定理6.1,6.5,6.7及定理6.8可分别推得推论6.3若线性多步法(22)满足根条件及条件P=0,。=。，则该方法是收敛的。推论64若线性k步法(22)满足根条件及(20),且Cp+i#0,则该方法是P阶收敛的。由于Runge-Kutta方法属单步法，则必零稳定。据此，推论6.2及定理6.3,6.7,6.8可分别推得S推论6.3若Runge-Kutta方法(29)满足%=1 则该方法收敛。，=i推论6.4若Runge-Kutta方法(29)满 足 阶 简 化 条 件。),。(),且夕为满足6 初及不等式p W m i n +1,2 +2的最大正整数,则该方法为p阶收敛的。6.4数值方法的有效实现一个微分方程数值方法构造出来后，要想真正在计算机上有效实现，求出合乎原问题要求的数值解，还需克服许多困难，如误差如何估计，步长及计算初值如何选取，隐式方法及多步方法如何处理等。在这些问题中，方法的阶与步长的选择是数值求解微分方程的首要问题，它直接影响计算精度和计算量。对于方法的阶，我们一般要求其不超过原问题真解的可微次数。但是，我们也必须注意到：高阶方法的稳定性能通常较差，因此阶数很高的方法一般不适用。对于方法步长的选取，从减少截断误差的角度来看，应采用小步长，但是步长取得过小，计算量将增大，而引起大的舍入误差。为此在步长选择时，我们必须二者兼顾，合理选取。通常可利用误差主项来确定步长，即要求所取步长使方法的误差主项不超过预定精度且接近5,据此在计算中调整步长。与此同时，我们也要求步长符合稳定性要求。数值方法主要分为显式方法与隐式方法，从计算实现来说，显式方法易于处理。为此，本节主要介绍隐式方法的有效实现，其主要倚重二类技巧，一是迭代技巧，二是预估-校正技巧。以下，我们将应用这二种技巧于隐或线性多步和隐或Runge-Kutta方法，仄而实现隐式行法有效求解常微分方程初值问题。6.4.1迭代技巧隐式方法的实现问题实质上是非线性方程的求解问题，因此第五章所介绍的迭代方法一般均可用于隐式方法的有效实现。本节以 Newton迭代法为例来说明隐式方法的实现。对 于 隐 式 线 性 多 步 法(2 2),不妨置 Qk=1，并 记 必=E (h 0 i fn+i-Q的+力，则该方法可写为2=0V n+k =+B k *0.(55)应用Newton迭代法于(55)有姆=此 k-卜 ia叫1 y%h“(tn+k,心-4,r=0,1,2,(56其计算终止准则可取为照 隽)-期/或帆空i/V陶+M腐)-=107-12)&err2 =107-12)r=y20 h*beta(3)*f(t2,y20)w;y21=y20 in v(eye(d)h*beta(3)*d f(t2,y20)*r;err 1 =norm(y21 y20);err2=norm(r);y20=y21;endy2=y20;tO=tO+h;yO=yl;yl=y2;end例6.4设有初值问题“y ,5 ，Wi(川V。同、7(0)=2，=0,(57)其中/，q 为已知函数，其使得该初值问题有唯一解x(i)=(2T T t)cos，y(t)=(2?r t)sin.).试取步长h=0.0 1,应用三阶Adams-Moulton方法(25)计算该初值问题，并作出其数值解的曲线图及整体误差图。解 应用算法6.2可获得初值问题(57)的数值解，其整体误差err=|Xn-X(in)|2及解曲线见图6.3,这里XR：=xn,X(tn)：=2 ,一(5)，tn=Tbh.由其整体误差图可知，该 Adams-Moulton方法关于初值问题(57)在求解区间0,5 内的计算是有效的，其精度达IO-量级。x 1 o7图 6.3 Adams-Moulton方法解问题(57)的整体误差图及解曲线图.对于隐式Runge-Kutta方法(2 9),其实现的关键是每步需解其sS个内部级值匕,7 2 =枷+八 aijf(tn+Cjh,Y n)(2=1,2,s)o 为j=i此，我们将其写成一个紧凑的形式YR=e 0 yn+h(A 01d)F(tn,Y n),(58)其中e=(M._ J)r,%=(略 工)，K i)6F(tn,Yn)=(f(tn+dh,Y i,n)T,f(tn+c2h,Y nf ,f(tn+c5h,YST)T.若应用N ew ton迭代法到(58),则得%+1)=力/)_ 卜以4 0 X严与产 t 设n e 0 yn h(A 0 Id)F(tn,Yr,=0,12 ，I其每步迭代初值Y f)可由同阶显式R un ge-K utta方法给出，计算终止准则可取为|K jr+1)-r)|或e 物6(4 =107-12)&err2 =107-12)F=f(t(l),Y 0(l);Y 0(2);f(t(2),Y 0(3);Y 0(4);r=Y 0 kron(on es(s,1),y O)h*kron(A,ey e(d)*F;d F =d f(t(l),Y O(1);Y O(2)zeros(2,2);zeros(2,2)d f(t(2),Y 0(3);Y 0(4Y l=Y 0i n v(ey e(s*d)h*kron(A,ey e(d)*d F)*r;errl =n orm(Y l Y O);err2=n orm(r);Y O=Y 1;en dy l=y O+h*kron(B,ey e(d)*L f(t(1),Y 0(l);Y 0(2)J);f(t(2),Y 0(3);Y 0(4)J)J;tO=tO+h;y O=y 1;en d例6 5长h=0.0 1,应用二级Gauss方法计算初值问题(5 7),作出其数值解的曲线图及整体误差图。解 应 用 算 法 6.3 可获得初值问题(57)的数值解，其整体误差err:=|Xn-X(tn)|2及解曲线见图6.4,这里Xn：=xn,7/n X(17 2)：=1(1n),g(九),tn=nh.由其整体误差图可知，二级Gauss方法关于初值问题(57)在求解区间0,5 内的计算是有效的，其精度达10-1量级。图 6.4 二级Gauss方法解问题(57)的整体误差图及解曲线图.6.4.2预估-校正技巧隐式方法的另一个有效处理途径是预估-校正方法。理论分析表明，迭代法中的初始迭代值若以与迭代公式同阶或低一阶的公式计算，则至多只需一次迭代，其迭代值精度即可与迭代公式本身的固有精度同阶。为此，我们可用与迭代公式同阶的显式公式先计算迭代初值(并称此显式公式为预估公式)，然后用迭代公式(称之为校正式)修正一次，即可获得其数值解。此方法称为预估校正方法。如我们以二步三阶显式公式(24)作为预估式，而以三阶Adams-Moulton公式(25)作为校正式，则司 得预校算法f预估-嬴+2=4 蜘+i+byn+2/z2%+i+.(60)校正：n+2=Vn+1+制5/陶+2,为+2)+8/n+l /J-其中启动值约以分别由原问题初值及同阶单步方法获得。为进一步提高预校算法的精度，我们通常以预估式及校正式的误差估计量来分别修正预估值和校正值。如在式(60)中，若取物+i x y(tn+1),yn x g(力n),则根据(21)有圾+2)-Vn+2%陶)，M+2)-Vn+2-)礴 (%)由上两式可得其事后误差估计式 4 1(力n+2)3 n+2 工(姗+2 为+2),沙陶+2)V n+2 一 工(纳1+2 一/+2),5 5(61)记/瓦+2,加计2,金+2分别为第+1个计算步的预估值，修正值及校正值，且分别以式(61)的右端作为预估值与校正值的修正量，则可得如下预估-改进-校正方法预 估：P n+2 =-4g九+i+oyn+2九(2%+1+fn),进改+&4-5Fn+2 f (tn+2,1m n+2),校 正:Cn+2=yn+l+卷(5&+2+8/n+l _/n)5改进：V n+2=C n+2 5(。计2 P rz+2),fn+2=/(力 九+2,枷+2).(62)上述方案中，计算预估式的改进值加计2时采用的是前一步的修正量累。计1-P计1)，其原因是当前步的校正值为+2此时还未算出。此外，在启动方案(62)时，我们可取。1-0 =0。算法6.4 预估-改进-校正算法function z=pmcm(a,b,yO,d,h)tO=a;cl=0 0 5;pl=0 0A=0 0 0;l/3 0 0;0 2/3 0;B=l/4 0 3/4;c=0 1/3 2/3 5;tc 1 =tO+c(1)*h;tc2=t0+c(2)*h;tc3=tO+c(3)*h;Yl=y0;Y2=yO+h*kron(A(2,1),eye(d)*f(tc 1 ,Y1);Y3=yO+h*kron(A(3,1:2),eye(d)*f(tc 1 ,Y1);f(tc2,Y2);yl=yO+h*kron(B(l:3),eye(d)*f(tc 1 ,Y1);f(tc2,Y2);f(tc3,Y3);for n=a+2*h:h:b11 =tO+h;t2=t0+2*h;p2=(-4)*y1+5*yO+2*h*(2*f(tl,y l)+f(tO,yO);m 2=p2+4*(cl-pl)/5;F2=f(t2,m 2);c2=yl+h*(5*F2+8*f(tl,y l)f(tO,y 0)/l2;y2=c2-(c2-p2)/5;f2=f(t2,y2);tO=t1;yO=yl;yl=y2;c 1=c2;pl=p2;end例 6.6取步长h=0.01,应用预估-改进-校正方法(62)计算二阶初值问题(2txt)+2x(t)=t3 In t,t e 1,5,l=1,V=0,(6 3)试作出其数值解的曲线图，并 求 出 其 全 局 误 差 max|功-OWngNX（tn）,这里 Xn N=400o解 记沙=，则初值问题(63)可化为(V(力)=y(t),t e 1,5,以力)=加力)一初+1ln。t G 1,5,(64)力=1,沙=0,其精确解为/=兴?2 In力-3)+7,y(t)=|t2(ln-7)+7.取步长h=0.01,应用算法(62)于初值问题(64)可得其数值解，该数值解外的曲线图见图6.5,其全局误差为err:=max xn x(tn)=6.7259e 008.0nN当用预校算法(60)计算初值问题(63)时，其数值解出的全局误差为err:max xn x(tn)=1.3945e 006.0nN此时，该算法精度低于相应的预估-改进-校正方法(62)的精度。16图 6.5 预估-改进-校正方法(62)应用于初值问题(63)的数值解.First Prev NeM/Last Go Back6.4.3刚性问题的数值处理考虑如下由Lambert22给出的常微分方程初值问题:问题I：1 A （2/i W Y f 2si n力-2 J y2（t）J1 2（cos力-si n力）J问题II：问题I 与问题II有同样的精确解t G 0,10（黝=（孀）=2ex p（T）（：）+（）V。,明表 6.8 显式Euler法求解问题I,I I 的整体误差利用显式Euler法并取步长h=0.1/2,(i=0,1,2,3,4,5,6)分别计算问题L II在末点tf=10处的数值解%(10)0 =1,2,3,4),其整体误差 j=|%(10)防(10)|(J =1,2,3,4)见表 6.8.h?：问题I问题 II功 300.18.7014e-0031.9646e-0021.7714e+1921.7679e+1950.054.1652e-0039.7393e-0030.0252.0377e-0034.8492e-0030.01251.0078e-0032.4196e-0030.006255.0118e-0041.2085e-0030.0031252.499 le-0046.0395e-0040.00156251.2478e-0043.0190e-0041.1388e-0041.1365e-004表 6.8表明：显式 Euler法求解问题I 是有效的；但该方法用于求解问题n 时，若取步长九=0.1,则产生巨大误差，若取步长h=0.05/7(i=0,1,2,3,4),则产生溢出，仅在步长逐步减半至非常 小 步 长 九=0.0015625时，算法才产生有效数值结果。对于上述数值实验，我们自然要问：同一数值方法求解形式上类似、且具相同精确解的二个初值问题时，其数值结果为何产生如此大的差异？显然从方法的相容阶角度我们无法解释上述现象，而只能从方法的稳定性来诠释其差异。为此，我们引入文献 21 中的一个稳定性结果.定理6.9设方法(44)的稳定域为S,d维线性系统|喘一四+G(65)I V (0)=甯的系数矩阵4 具 特 征 值 方 法(44)求解具不同初值的系统(65)所得数值解为 枷,许,则lim(yn-zn)=0(66)九 一oo当且仅当hXi e S,i =1,2,(67)在条件(67)下，(66)式表明数值误差在传播过程中保持有界,且随n无限增大而趋于零，按此意义我们视计算过程是数值稳定的。由定理6.9,若要显式E ul er法分别关于问题L I I 数值稳定的,则依次应有九MGSEE=N GC:|N+1|1 ,i=1,2;(68)hX1 e SE E=zeC:z+l 1,2 =1,2,(69)这里，SE E为显式Euler法的稳定域，而X=-1,夭2=-3,Af=-1,对=-1000(70)为问题I,I I 的系数矩阵的特征值。将(70)分别代入(68),(69)可解得20 h -0.667(关于问题 I)(71)O0 h P=卷 匐 （号）|阳 因 6凡（誓其中p（），%（）分别为矩阵的谱半径与特征值。因此，对于刚性问题”3）,其 Lipschitz常数满足L Lmn=sup|吗河 2 sup 阴吃（吗y）|1.n）eD （t,y）eD 口少 7由此可见,非线性刚性问题具有异常巨大的经典Lipschitz常数的特征。人们可以断言：刚性是初值问题本身所固有的，其实质是数值求解慢变解时，存在快速衰减的扰动，这种