线性多步法续

1、第4章线性多步法4.1线性多步法的一般公式前面给出了求解初值问题(1.2.1)的单步法，其特点是计算 加 时只用到八 的值，止匕时不。，乃,，八J居 的值均已算出.如果在计算乂+1时除用人 的值外， 还用到的值，这就是多步法.若记取=/+粉,h为步长， 1n巾则线性多步法可表示为i-12:，(4.1.1)i-0i-l其中名，以为常数，若以3+f3Ho(即不同时为零)，称(4.1.1)为线性k步 法.计算时用到前面已算出的k个值%，儿】，居-5.当/“0时，(4.1.1)为显 式多步方法，当则称(4.1.1)为隐式多步法.隐式方法与梯形方法一样，计 算时要用迭代法求M+1.多步法(4.1.1)的

2、局部截断误差定义也与单步法类似.举例来说，对于初值问题y' = -y+x+1, y(0)=1,步数k=2时，线性多步法 表小为yn 1 = ： 0yn： Unh(：fn 1：0力：$)，n =1,2,当员1 0时，格式为显示的：yn 书=a0yn +%yn+hBo(yn +Xn+1)+B1(yn/+Xn+1, n = 1,2,，而凤】H°时，格式为隐式的：%1 =:0%：1%1h只(-1Xn11)：。(-%Xn1)1(-%/ %二 1, n =1,2,。定义4.1设y(x)是初值问题(1.2.1)的精确解，线性多步法(4.1.1)在外+1处的局部截断误差定义为jt-ibi-一

3、 - 上一，'二二.(4.1.2)若丁/1=°5"),则称线性多步法(4.1.1)是p阶的.如果我们希望得到的多步法是 p阶的，则可利用Taylor公式展开，将或+】在 。处展开到状41阶，它可表示为“qm +CM&)+ C33”/)+ %也兴+】门叫/)+8炉力（4.1.3）注意，（4.1.2）式按Taylor展开可得月-1hTa+i = M/ +以淡）血工用y区3=0i=-l=） + 如'（/） + 尸”'（/） + （:十罩1/*）（/）+ MTfj我）-友£ wx）一期1（/）+-g-“（0）+,】经整理比较系数可得Cq =

4、，二Q、G£(t)% -工乩 r=Dis-1 I Jt-1C厂不1一豆(7)%一小豆5尸友八空,“+1若线性多步法(4.1.1)为p阶，则可令=孰=,=。/。,4+产0于是得局部截断误差-(4.1.5)右端第一项称为局部截断误差主项.C>+1称为误差常数.要使多步法(4.1.1) 逼近初值问题(1.2.1),方法的阶p> 1,当p=1时，则Co=C=O,由(4.1.4)得Jfc-l Jt-L/十为十+ %i T，E肛一2自二T（4.1.6）称为相容性条件.公式(4.1.1)当k=1时即为单步法，若内-1 = 0,由(4.1.6)则得% = L d=1式（4.1.1）就是+

5、%，即为Euler法.此时。口 =；=0,方法为p=1阶.若足户0,由。口 = °得 =1,为确定尸4及凡，必须令G=q = 0,由（4.1.4） 得1.1 =1及/;£此时（4.i.i）就是八八*a 十九）即为梯形法. 由 一 一 ，:一 ' 一一 ' _:31J t> 41212故p=2,方法是二阶的，与3.1节中给出的结果相同.实际上，当k给定后，则可利用（4.1.4）求出公式（4.1.1）中的系数%及以，并求得Tx+i的表达式（4.1.5）.1.2 Adams显式与隐式方法形如jt-i一 丁， , '' 以1一、（4.2.1）的

6、k步法称为Adams方法，当#-L。时为Adams显式方法，当员】工。时, 称为Adams隐式方法.对初值问题（1.2.1）的方程两端从4到/u积分得双/+1）闪/）=飞九式初办显然只要对右端的积分用插值求积公式，求积节点取为或衿国，2即可推出形如（4.2.1）的多步法，但这里我们仍采用Taylor展开的方法直接确定 （4.2.1）的系数弁G =-L0,代-1）.对比（4.1.1）可知，止匕时%-0,只要确定员1，尸口，尸即可.现在若k=4且乩I =0, 即为4步的Adams显式方法八* =居+见优人+回兀.1十匹手一十自£与）其中乩,尸卜为，区为待定参数，若直接用(4.1.4),可

7、知此时ci-% =1-1- o自然成立，再令c1 = cq =c4 =o可得区+回+用+饱=1 022-3仇=-；iui +4氏+9屎 从7凡+27&=-解此方程组得用=25卜福心胃心弋.由此得到251720=乂 +与3以79以+37九-9九)(4.2.2)(4.2.3)于是得到四阶Adams显式方法及其余项为Ml北=卷炉ySG斗口6。)r乙"若反。，则可得到p=4的Adams隐式公式，则k=3并令GGqc;o由(4.1.4)可得/_】+ + £i+同=1 1/_1+ 优 +% = JBH瓦J1解得员一方网嗜血一盘血 W，而一起5/藤, 于是得到四阶Adams隐式方

8、法及余项为：二." 二，- -.一：.(4.2.4)丁+|-；二'；:(4.2.5)r riu般情形，k步Adams显式方法是k阶的，k=1即为Euler法，k=2为= y戛工-<k=3 时,丁科1 二 K+'（23/,-16 九+5冗k步隐式方法是（k+1）阶公式，k=1为梯形法，k=2为三阶隐式Adam公式 乂 + （5/1+1 + 也-Xh1）k步的AdamspJ法计算时必须先用其他方法求出前面 k个初值Sq才能按给定公式算出后面各点的值，它每步只需计算一个新的f值，计算量少，但 改变步长时前面的-卜乂/、居心也要跟着重算，不如单步法简便.例4.1用四阶显

9、式Adams方法及四阶隐式Adams方法解初值问题产_尸+工十1”式41,“0）=1,步长h=0.i用到的初始值由精确解尸（区）=仁”+再计算得到.解 本题直接由公式（4.2.2）及（4.2.4）计算得到.对于显式方法，将汽孔尸）=-，+彳+1直接代入式（4.2.2）得到九+居+石（55工79人.】+ 37力_9人.）”3.4,.9其中','1对于隐式方法，由式（4.2.4）可得到=Z1 +学网”十1】+ D+1" - 5£_】+九直接求出八*1,而不用迭代，得到Q1%+1 =六十赤父9“+1 + 9 + 19为一5人“十人一力用=2,3,，90-324 y

10、计算结果如表所示.表4-1 Adams方法和Adams隐式方法的数值解与精确解比较精确解贝题尸£,+%Adams法Adams隐式方法%网而>珈1先0.3L040S1S221040818012.lxl(rTD.41.070320051 07032252187X10-61 070319663 9x107151.106530661.106535484 82x1(1.106530145 2x1070.61.148811641 148818416.77x10-61448811016.3x10?071.193585301.196593408.10x10*1496504597.1x1Q.3L

11、2493羽961.249338169.20x101 249328197.7X1O10.91.306569661.306579629.96X1O(5306568848.2X1O1101.367879441 367889961 Ui M1,367878593.5x1a11.3 Adams预测-校正方法上述给出的Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由 于每步要做迭代，计算不方便.为了避免迭代，通常可将同阶的显式 Adams方法 与隐式Adamsf法结合，组成预测-校正方法.以四阶方法为例，可用显式方法 (4.2.2)计算初始近似乂2,这个步骤称为预测(Predictor),以P

12、表示，接着计算f值(Evaluation),这个步骤用E表示，然后用隐式公式(4.2.4)计算入+1,称为校正(Corrector),以C表示，最后再计算)*1 =/(4“,八+1),为下一步计算做准备.整个算法如下：血预测；p： y=八+五3M -59九】+ 37fz-对小(4.3.1)求值* Et*=(演+1，y3)校正c： X " +5。益19力-5九十九)求值比%1 = /(%1必+1)公式(4.3.1)称为四阶Adams预测-校正方法(PECE).其matlab程序如下function y = DEYCJZ_adms(f, h,a,b,y0,varvec,type) for

13、mat long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(i) = y0;y(2) = y0+h*Funval(f,varvec,x(1) y(1);y(3) = y(2)+h*Funval(f,varvec,x(2) y(2);y(4) = y(3)+h*Funval(f,varvec,x(3) y(3);for i=5:N+1v1 = Funval(f,varvec,x(i-4) y(i-4);v2 = Funval(f,varvec,x(i-3) y(i-3);v3 = Funval(f,varvec,x(i-2) y(i-2);v4 = Fun

14、val(f,varvec,x(i-1) y(i-1);t = y(i-1) + h*(55*v4 - 59*v3 + 37*v2 - 9*v1)/24;ft = Funval(f,varvec,x(i) t);y(i) = y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24;endformat short;问题1(4.3.1 )中的第四步在程序中那一行实现了?利用(4.2.2)和(4.2.4)的局部截断误差(4.2.3)和(4.2.5)可对预测-校正方 法(4.3.1)进行修改，在(4.3.1)中的步骤P有yd/ - y：中之言个丫720问题2这个约等式是怎么得来的？什么方法？对

15、于步骤C有P 195 (5)yn i - yn 1 h y (%)720两式相减可得于是有U5(/5 )、720 / p、h y (Xn) : "270 "n 1 一 yn 1)251yLi) * -旃 33 - kG19yRG -yi % 痂必i -温若用/工代替上式八+】，并令19270(端-端)显然K1为川比此1，式+】更好,但注意到尸言的表达式中吃】'是未知的，因此改卜面给出修正的预测-校正格式(PMECME).P :匕1 ="十元(55人59人37力(4.3.2)股溜二匕+券w-2C：其Ta店+19力-5力九)2419财'- ym+l J

16、ki .7口 Ok+i -齐+i) u ¥ U1产：力+1 =y(Xjt+lrL*+l)经过修正后的PMECME式比原来PEC弗式提高一阶.问题3试着编出该程序！ ！4.4 Milne 方法与Hamming方法与Adams显式方法不同的另一类四阶显式方法的计算公式形如（4.4.1）这里/o,4黑口为待定常数，此公式也是k=4步方法，即计算八+1时要用到居,tMtJt，4个值.为了确定同，邑L当然可以利用公式（4.2.1）直接算出，但下面我们直接利用Taylor展开式确定，凤，自，伉，使它的阶尽量高.方法（4.4.（1） 部截断误差为T"i =尸（/十川+双/一+从昆尸）十乩

17、尸乂/一砥十足力工/一展）将它在五点展成Taylor级数，得% = y（三）+-A 铝 E+纱 O 勺尸阳+2产+-«）-湖气）+ 萼y“（x*）-萼V"（%） + 粤/）&）-粤*（“）+ ,】-（h2h31%* （&） +月LGJ -匕。*）+彳。'（4）-臼。'（/”丁，（/）-】+日I "、n? iiz y （2心（2以）M （2血）* I1产吗）-23口苛*（/）+彳卜L=1 + 3- 仇十4+）历'1） + g g+同+ 2油3小） +12714 口、下川/ 、,1*118 a j 4 （41 /、（7+不用5户

18、。用尸（/）+一力+ W4+ £/。诙¥（勺）十12510*可一祕一下耳）.叫/）+。（盾要使公式的阶尽量高，要令前 3项系数为0.即2 g4一（乩十回十2）=。4十月+2乩=0,四+4=3848i* V（*?/T 解得;2=亍/1=-§，凤=5代入公式， 3的系数为0,故,一；:'.：际'（4.4.2）于是得四阶方法一4-!|1-:、1（4.4.3）称为Milne公式，它的局部截断误差为(4.4.2).与(4.4.3)配对的隐式方法为k=3的多步法，它的一般形式可表示为"八 + iiJS" + 牝乂c + M0-1力*+ 员&

19、lt; + A/B-1)要求公式的阶p=4,可直接用(4.2.1),并令7 = Q = j = q = 0,可得IQ +% + 劭 T-曲- 2牝+氏+屈+凡=1d 4 + 4% + 2a队=(4.4.4)-勺-犯+ 3储+3自=1 鼻1 + 16% +4/_i -4/?i - 1 L14若令心 =°,可解出& = 0q=，/"=四=才风=可，于是得到下列四阶方 法.用:”.1 .；,-i(4.4.5)称为Simpson公式，它的局部截断误差为-一二(4.4.6)用Simpson公式与Milne公式(4.4.3)相匹配，用(4.4.3)做预测，(4.4.5) 做校正

20、，由于(4.4.5)的稳定性较差，因此通常较少使用.为了改善稳定性，可重新选择四阶的隐式公式，Hamming!过试验，发现在(4.4.4)中若令% = °,得到_ 9_ 1_3 _ _6 _3的公式稳定性较好，此时(4.4.4)的解为1。=§&=一+*="尻=*仄=一豆 , 于是得四阶多步法13，.:,一 .i 1；< 1(4.4.7)称为Hamming式，它的局部截断误差为工. 上，工厂(4.4.8)用Milne公式(4.4.3)与Hammin也式(4.4.7)相匹配，并利用截断误差公式 (4.4.2)与(4.4.8)改进计算结果.4-4yn 1

21、= ynq ;(2 fn - fn- 2心)3(4.4.7，)1yn 1 =819yn fl 3h(fn1 2fn该算法称为Hamming®测-校正法。类似Adam颉测-校正格式（4.3.2）,可得以 下的修正的 milne-Hamming预测-校正格式：4产尸幻=Aq +可双2式一工5+2二）119(4.4.9)m y黑次乙】+行-统f）百乜警=（小,置）13C A；+i =石9八-y）+ 坪点 + 2工一人.j OO9M;以荷=下工1 工i 一3工1）产,/T+1 j（Xji+1，Am+】）附hamming程序function y = DEYCJZ_hm(f, h,a,b,y0,

22、varvec,type) format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,x(1) y(1);y(3) = y(2)+h*Funval(f,varvec,x(2) y(2);y(4) = y(3)+h*Funval(f,varvec,x(3) y(3);if type = 1 %hamming 预测校正法for i=5:N+1v1 = Funval(f,varvec,x(i-3) y(i-3);v2 = Funval(f,varvec,x(i-2) y(i-2);v

23、3 = Funval(f,varvec,x(i-1) y(i-1);t = y(i-4) + 4*h*(2*v3 - v2 + 2*v1)/3;ft = Funval(f,varvec,x(i) t);y(i) = (9*y(i-1) -y(i-3) +3*h*(2*v3 + ft-v2)/8;endelse % 修正的hamming预测校正法p0 = 0;c = 0;for i = 5:N+1v1 = Funval(f,varvec,x(i-3) y(i-3);v2 = Funval(f,varvec,x(i-2) y(i-2);v3 = Funval(f,varvec,x(i-1) y(i

24、-1);p = y(i-4)+4*h*(2*v3 - v2 + 2*v1)/3;M = p - 112*(p0 - c)/121;F = Funval(f , varvec, x(i) ,M);c = (9*y(i-1) -y(i-3) +3*h*(2*v3 + F-v2)/8;y(i) = c + 9*( p - c)/121;p0 = p;endendformat short;例4.2用四步四阶显式Milne公式及三步四阶隐式Hammin也式解初值问题l一 + m + LQ & x S 巾1步长 h=0.i初值弘乃仍由精确解迎）=尸+工给出，要求计算到J 二03为止，给出计 算结果

25、及误差，并与例4.1结果比较.解直接用公式（4.4.3）及（4.4.7）计算.用Milne法计算公式为04%+k + x QA -三 3,4其中， ,.1 = 一%十加十1 二 0工-71 + + 1 = -1 00483742 + 1 1 - 0.09516258f2 =$+旃 +1 = -1.01873075+ 1 2 = 0.18126925£ 工-券 +4+1 = -1 04081822 + 1.3 = 0,2591317804筋=R + 亍 x （2力-力 + 2力）=1 07032260八三一M +x4 +1 = 0.3296774004-三力=M + 乂（2从-/ +

26、2力）=1 10653229误差I y（x，-y41= 2.55x1V,卜（工%）一梅卜 L63*10f用Hamming法（4.4.7）计算公式为%+1 =9。笫-%十?&九+ 2M -九）103=gX （沙国一 Al） +目*仁八斗冬川+1斗2/, 一人_）可解得，n=2,3,4y.+1 = x一八7） + x（+i +2/, 一/7 o. Jo. J103,为= 7*（9凡-（电+1 + 2力-西）0.36. J-X 8.168576 + x r 56737592 = 1 040818028.38.3八=一内十西-1 = 0.2591819810 3M =白乂（9%-冯）+三 （演

27、7 + 2人-入户1 07031966八-乂 +。+1 =0.329680343 - Mg仍）+ 4 乂5 + 】+ 24 - 1.10653010误差卜（曲）= 2,0xl0y4） -y4| = 3,9xl0-7|y（x5） -y5| = 5.6x10从所得结果可见Milne方法误差比显式AdamspJ法误差H&小，而Hamming 方法与隐式Adams方法误差相当.例4.3将例4.2的初值问题用修正的Milne-Hamming预测-校正公式计算5及汽,初值冲，乃，乃，乃仍用已算出的精确解，即%=5 = 1 00483742r72 = L0187307工1=13081822给出计算结

28、果及误差.解 根据修正的Milne-Hamming预测-校正公式（4.4.9）得y； '1 10653299V严=W 十皆乂；-X）-1.106530364A.3/）*' + 均+ 1 = 0.393469636IIn 7.-9+?川力“尹十"力】 OO-1.106530419冷=另一言另一火）=1 10653061力=-尸5 +#5 +1 = 0,3946939M =乃+子气2工-£ + 26=1.14881135靖=一产产十 % +1 = 0 45118865In 3y： = - &9一加 + = x+ 2工元）=148S114488n =X 一

29、言*（X -yl） = l 14密ii6iI一匕 1= 497 父1L3g） -4 |= 2.61x10从结果看，此方法误差比四阶Adams急式法和四阶Hammin/法小，这与理 论分析一致.讲解：线性多步法（4.1.1 ）的局部截断误差定义为与单步法相似，可表示为（4.1.2 ）,即fc-1Jt-1小 r氏+m二叩氏访）tZ i=0i=-l只要直接将右端各项在 2处展成Taylor公式，根据公式阶数为P阶，即北川二0（加附）按心的幕整理，令/次1必立各项系数为0,则可求得相应的线性多步法 及其局部截断误差，这里只用到一元函数的Taylor展开.因此不必记系数满足的 公式（7.5.4 ）,只要

30、直接展开即可，它不但可以求出 Adams显式与隐式公式以 及Milne公式，Hamming式等，还可以求出任何需要的多步法公式，下面再给 出两个例题，说明如何直接用 Taylor展开的方法.例4.4 解初值问题v'=fMy）,以刈）二的用显式二步法=4% +%7加1+双用工+ A/m-i）,其中4二工!匚/（/-卜居-1）.试确定参数% 以卜用使方法除数尽可提高.并求局部截断误差.解 本题仍根据截断误差定义，用Taylor展开确定参数满足的方程，由于-=4（% +A）-3（/）一叩（% -A）-型如+如1/-A） 京族3二量+加气）+, + y''' （/） + 了伙/）+。面）-%y- 现权（Q -印3g）-即+gv区）-*收）十。册）1 1 .二（1 一 % %»（/）+。+ / _ 凤%的气）+ - + 龄卜勺（） + 伊34必,"' g+（5 -噌京 必产+。的为求参数/，” 自，4使就地介数尽量高，可令1 一4一/=0, 1+的一昆一/1=0 ；一；%+d=0, ! + !。；