常微分方程数值解

1、第四章常微分方程数值解课时安排6学时教学课型理论课教学目的和要求了解常微分方程初值问题数值解法的一些基本概念，如单步法和多步法，显式和隐式，方法的阶数，整体截断误差和局部截断误差的区别和关系等；掌握一阶常微分方程初值问题的一些常用的数值计算方法，例如欧拉(Euler)方法、改进的欧拉方法、龙贝-库塔(Runge-Kutta)方法、阿达姆斯(Adams)方法等，要注意各方法的特点及有关的理论分析；掌握构造常微分方程数值解的数值积分的构造方法和泰勒展开的构造方法的基本思想，并能具体应用它们导出一些常用的数值计算公式及评估截断误差;熟练掌握龙格-库塔(R-K)方法的基本思想，公式的推导，R-K公式中

2、系数的确定，特别是能应用“标准四阶R-K公式”解题；掌握数值方法的收敛性和稳定性的概念，并能确定给定方法的绝对稳定性区域。教学重点与难点重点：欧拉方法，改进的欧拉方法，龙贝-库塔方法。难点：RK方法，预估-校正公式。教学内容与过程4.1 引言本章讨论常微分方程初值问题(4.1.1)l$a)=几的数值解法，这也是科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法.通常我们假定(4.1.1)中f(x,y)对y满足Lipschitz条件，即存在常数L0,使对1必已应，有(4.1.2)于SJ一了(兀旳)卜庞一旳则初值问题(4.1.1)

3、的解存在唯一.假定(4.1.1)的精确解为畑,求它的数值解就是要在区间禺切上的一组离散点=上求用)的近似儿仍，必-.通常取叫=&+曲伍=0丄)，h称为步长，求(4.1.1)的数值解是按节点确=12)的顺序逐步推进求得旳：首先，要对方程做离散逼近，求出数值解的公式，再研究公式的局部截断误差，计算稳定性以及数值解的收敛性与整体误差等问题.4.2 简单的单步法及基本概念4.2.1 Euler法、后退Euler法与梯形法求初值问题（4.1.1）的一种最简单方法是将节点的导数讥）用差商y（心+?-y（心）代替，于是（4.1.1）的方程可近似写成M召+i）珀y（xn）+hf（xn,y）（xri）,n=0?

4、1?-从坯出发$2）7显）=Vo,由（4.2.1）求得沪%+防U,丹）=戸再将戸-巩皿代入？.1）右端，得到讥花）的近似/厂戸+紅仏乃），一般写成(4.2.2)称为解初值问题的Euler法.Euler法的几何意义如图4-1所示.初值问题（4.1.1）的解曲线y=y（x）过点%（忑0），从坯出发，以畑肌）为斜率作一段直线，与直线“心交点于珂显然有戸=兀+幼坯,丹），再从珥出发，以/仪11）为斜率作直线推进到兀=心上一点3,其余类推，这样得到解曲线的一条近似曲线，它就是折线占乙Euler法也可利用讥心+1）的Taylor展开式得到，由y(n+h)=y(Aj+hylUJ+#l(OE(心耳+1)(4.

5、2.3)略去余项，以儿讯耳），就得到近似计算公式（4.2.2）.40另外，还可对（4.1.1）的方程两端由入到+1积分得m耳+i）一儿=b+i/（兀巩X）必(4.2.4)若右端积分用左矩形公式，用儿兀），儿+1（+1）,则得（4.2.2）.如果在（4.2.4）的积分中用右矩形公式，则得M+i=入+紅（禺+1，尹肝1），挖=0，(4.2.5)称为后退（隐式）Euler法若在（4.2.4）的积分中用梯形公式，则得人+1=人+了氏+1，人+1儿找=丄（4.2.6）称为梯形方法.上述三个公式(4.2.2),(4.2.5)及(4.2.6)都是由町十算儿+】，这种只用前一步即可算出儿+1的公式称为单步法，

6、其中(4.2.2)可由兀逐次求出丹丿2,的值，称为显式方法，而(4.2.5)及(4.2.6)右端含有了(兀小儿+1)当f对y非线性时它不能直接求出儿+1,此时应把它看作一个方程，求解叫+1，这类方法称为稳式方法此时可将(4.2.5)或(4.2.6)写成不动点形式的方程兔+i=加#(+1，儿+i)+g这里对式(4.2.5)有0=h盘,对(4.2.6)则居二儿+/(J，g与儿+1无关，可构造迭代法讯节=诃(耳+2鹽)+g牛0丄(4.2.7)由于煮兀刃对y满足条件(4.1.2),故有(5+1)充邛一儿+1卜啊乳耳+1，理J一了(心+1必+1)“加忖?1-儿+1|当hL或1-，迭代法(4.2.4)收敛

7、到儿+1,因此只要步长h足够小，就可保证迭代(4.2.4)收敛对后退Euler法(4.2.5),当h时迭代收敛，对梯形法(4.2.6),当曲彳时迭代序列收敛.例4.1用Euler法、隐式Euler法、梯形法解y1=-y+x+l?y(0)=1取h=0.1,计算到x=0.5,并与精确解比较.解本题可直接用给出公式计算.由于瓦凡巧=_y+x+l,h=O.l?xo=Ojq=1,Euler法的计算公式为儿+i=儿+h（yn+心+1）=（1-却儿+hxn+h=0.9兀+0.1+0.1n=0时，戸=叽+01000000其余n=l,2,3,4的计算结果见表4-1.对隐式Euler法，计算公式为+1=入+血（一

8、+1+疋話+1+1）解出儿+1=匕仇+叽+曲）=+以+1心+0.11）当n=0时，戸=占（必+0J可+011）=1-009091.其余n=l,2,3,4的计算结果见表4-1.表4-1例4.1的三种方法及精确解的计算结果為EulerStyn隐式Euler法yn梯形法精确解y（鬲）011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.120

9、9221.1062781.106531对梯形法，计算公式为人+1=儿+|（+1）+（一几+1+兀+i+1）解得当n=0时，戸=肓1.9+0.21)=1.004762,其余n=l,2,3,4的计算结果见表4-1.本题的精确解为F(瓦)瓦+巳，表4-1列出三种方法及精确解的计算结果.4.2.2单步法的局部截断误差解初值问题(4.1.1)的单步法可表示为入+1=y+h帆耳，人，人+1屈用=0丄(4.2.8)其中妙与有关，称为增量函数，当少含有儿+1时，是隐式单步法，如(4.2.5)及(4.2.6)均为隐式单步法，而当不含沪1时，则为显式单步法，它表示为(4.2.9)如Euler法(4.2.2),矶兀

10、y,h)仏y).为讨论方便，我们只对显式单步法(4.2.9)给出局部截断误差概念.定义2.1设y(x)是初值问题(4.1.1)的精确解，记(4.2.10)Tn+1=-尹(耳)h(xn,y(xn),h)称为显式单步法(4.2.9)在+】的局部截断误差.T+i之所以称为局部截断误差，可理解为用公式(4.2.9)计算时，前面各步都没有误差,即Vn=X)，只考虑由心计算到心+i这一步的误差，此时由(4.2.10)有Xn+l)-yn+l=心+1)一儿+血帆，儿周=巩耳+1)一巩心)一应就耳，巩耳)，閒=監+1局部截断误差(4.2.10)实际上是将精确解y(x)代入(4.2.9)产生的公式误差，利用Tay

11、lor展开式可得到丫曲=(血?).例如对Euler法(4.2.2)有帆和=/(),故瞪+1=肌耳+J-肌心）-紅（耳丿（）”（）+0（妒），称若显式单步法(4.2.9)的局部截断误差=畀（耳+h）-y（xn）-hy（xK）=押#（心）+押#11耳）+=0（小它表明Euler法(4.2.2)的局部截断误差为几+i为局部截断误差主项.P+1定义2.2设巩巧是初值问题(4.1.1)的精确解,),芒是展开式的最大整数，称F为单步法(4.2.9)的阶，含汕的项称为局部截断误差主项.根据定义，Euler法(4.2.2)中的戸=1故此方法为一阶方法.对隐式单步法(4.2.8)也可类似求其局部截断误差和阶，如

12、对后退Euler法(4.2.5)有局部截断误差T話+i=尹(心+1)一巩一幼(G+1J(叫+J)=只心+毎一巩耳)一劭(耳+曲)=切优)+?”优)+牛門耳)+-心1耳)+妙(耳)+故此方法的局部截断误差主项为(耳)，芒=1，也是一阶方法对梯形法(4.2.6)同样有T+i=讥3）-巩耳）-?/（心（耳）+/（耳+i耳+1）=小+却-巩心）-（耳+血）-护乜）+。国）它的局部误差主项为-巨卩（兀），戸=乙方法是二阶的.4.2.3改进Euler法上述三种简单的单步法中，梯形法（4.2.6）为二阶方法，且局部截断误差最小，但方法是隐式的，计算要用迭代法为避免迭代，可先用Euler法计算出+1的近似儿+

13、1,将（4.2.6）改为(4.2.11)(4.2.12)称为改进Euler法，它实际上是显式方法.即儿+1=儿+才了（耳必）+/K+1必+皆（耳必）右端已不含氐+1.可以证明乙+=0（妒），戸=2,故方法仍为二阶的，与梯形法一样，但用（4.2.11）计算入+1不用迭代.例4.2用改进Euler法求例4.1的初值问题并与Euler法和梯形法比较误差的大小.解将改进Euler法用于例4.1的计算公式人+i=儿+耳+1)+(一几_h(yn+1)+1+DI叩-宇反+忖兀紳”+“+)=0.905+0.095xh+0.1当n=0时，乃=0加恥+。风+0.1=1-005000.其余结果见表4-2.表4-2改

14、进Euler法及三种方法的误差比较耳改进E皿广法凡误差|y吐-y(xn)|Eukt方法-y(xn)|梯形法|yn-y(xn)|0.11.0050001.6xW44.8xl037.5xl050.21.0190252.9xl0-48.7xl031.4x100.31.0412184.0x101.2X10-21.9xl040.41.0708024.8xl0-41.4xl0-22.2x100.51.1.070765.5xW4l.SxW22.5x10从表4-2中看到改进Euler法的误差数量级与梯形法大致相同，而比Euler法小得多，它优于Euler法.讲解：求初值问题（4.1.1）的数值解就是在假定初值

15、问题解存在唯一的前提下在给定区间禺切上的一组离散点0=%“/y）=-y+Io=oo=i=oi于是当n=0时俎=了起,兀）=一坯+眄+1=0為=了（础+妮兀+争1）=（一1+*）必+1一）可+1=02居=畑+”屁+怠）=（1+曲牛+討+（1一扌+牛可+1=0.04右氐4=/（0+局兀+|可h2h3h1h3=（-l+-_+_o+（l-+-_）Xo+1=0.09525于是丹=兀+（組+生+%+屁）=1+*0,使对F忆已丘，均有I帆爲久-帆爲z?h)|L|y-z|则方法(4.3.2)收敛，且亡厂定理证明略.可见3.4.4.2 绝对稳定性用单步法(4.3.2)求数值解沪戸，小榔，由于原始数据及计算过程舍

16、入误差影响,实际得到的不是几而是久二儿*仇，其中久是误差，再计算下一步得到以Euler法为例，若3厂儿-儿，则(4.4.1)如果卩+肚y卜1,则从十算到九+i误差不增长，它是稳定的但如果条件不满足就不稳定.例4.4y=-100y,y(0)=1,精确解为畑=严，用Euler法求解得儿+1=儿+纱区必)=(1一10陶儿用丄若取h=0.025,则儿+1一匸讥，当心孔加儿l=co,而巴)吧严”=0,显然计算是不稳定的.如果用后退Euler法(4.2.5)解此例，仍取h=0.025,则片+1=入+VM+i=儿_2二片+i,即尹曲=丄儿显然当，计算是稳定的.由此看到稳定性与方法有关，也与几有关，在此例中几

17、=一1在研究方法的稳定性时，通常不必对一般的f(x,y)进行讨论，而只针对模型方程(4.4.2)y1=2y?Re(2)0这里兄可能为复数规定是因为0时微分方程(4.4.2)本身是不稳定的，而讨论数值方法(4.3.2)的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题(4.1.1),若将f(x,y)在氐血)处线性展开，可得廿儿)+耳，儿)(x-耳)+于；耳-儿)于是方程(4.1.1)可近似表示为它表明用模型方程(4.4.2)是合理的，至于模型方程(4.4.2)中所以用复数入是因为初值问题(4.1.1)如果是方程组，即ER贝伸是(皿乂皿)阶矩阵，其特征值可能是复数当然对单个方程，入

18、就是实数，此时只要规定几0即可.用单步法(4.3.2)解模型方程(4.4.2)可得到(4.4.3)其中依赖所选方法，如用Euler法则(4.4.4)入+i=儿+=1+泌)儿，運(血刃=1+血兄此时由(4.4.1)看到误差方程也为仇+1=锁泌)必，与(4.4.4)是一样的.因此对一般单步法(4.3.2)误差方程也与(4.4.3)一致.下面再考虑二阶R-K方法有对四阶R-K方法，可得定义4.2将单步法(4.3.2)用于解模型方程(4.4.2)，若得到(4.4.3)中的丨E兄)|1则称方法是绝对稳定的在复平面上复变量h几满足丨E他乳)|1的区域，称为方法(4.3.2)的绝对稳定域，它与实轴的交点称为

19、绝对稳定区间.例如对Euler法，|E(M|=|1+1U|1在复平面h几上是以(-1,0)为圆心，以1为半径的单位圆域内部，当几为实数时，则得绝对稳定区间为-2应乳0,因几0,故有在例4.4中h=0.02时方法稳定，而例中取h=0.025故不稳定.对后退Euler法(4.2.5)，儿+】=儿十力叽因几0方法都是绝对稳定的.二阶R-K方法的绝对稳定区间为-2.泌0.三阶R-K方法的绝对稳定区间为彳力四阶R-K方法的绝对稳定区间为_27S5例4.5用经典四阶R-K方法计算初值问题=-20X0x1)(0)=1步长取h=0.1及0.2,给出计算误差并分析其稳定性.解本题直接按R-K方法(4.3.12)

20、的公式计算因精确解为焜=厂血其计算误差比朋少如表所示.0.20.40.60.81.0h=0.10.0927950.0120100.0013660.0001520.000017h=0.24.9825.0125.0625.03125.0从计算结果看到，h=0.2时误差很大，这是由于在入=-20,h=0.2时入h=-4，而四阶R-K方法的绝对稳定区间为-2.485,0,故h=0.2时计算不稳定，误差很大.而h=0.1时h几=-2,其值在绝对稳定区间-2.485,0内，计算稳定，故结果是可靠的.讲解：由于微分方程初值问题数值解公式求出的是一个逐次递推的过程，因此原始数据误差及计算过程舍入误差对解的影响

21、就是数值方法绝对稳定性研究的问题，如果由十算7話+1误差不增长，方法就是绝对稳定的。为使问题得到简化通常就是将方法用于解模型方程(4.4.2)，对于单步法得到的差分方程为儿+1二運佻几)儿，由于模型方程的于(兀旳=勒，代入Euler法儿+=儿+虬二兀+嘉厲=(1+曲刃兀，得E(hZ)=+hl，对二阶r-k方法，例如，用改进Euler法儿+1=人+专1/孔必)+了耳+曲必+幼(兀,儿)=儿+轨+几(儿+曲恥)=1+泌+浑-儿应（泌）=1+泌+*3兄尸E（h2：=+hl+-（泌）2+-（加对三阶R-K方法有2印E(hZ)=+hl+-(M)2+-(閃+-(血兄对四阶R-K方法有23!41只要国血刃|

22、cl方法，就是绝对稳定的，这时的值当n增大式是减少的，故计算稳定。这时舍入误差影响可忽略不计，而当国谀刃卜】，则增大，方法不稳定，计算结果是不可靠的。因此用显式单步法必须使回认刃11,也就是步长选择要满足这一要求。对于隐式的梯形公式儿+i二儿+文LA，儿）+了（心小儿+1）.将模型方程*，即n）=代入得仏=儿+日恥+叽1+泌于是儿氣二_泌儿注意丸，于是有|1+扣|欧冈*n，对gg成立。这就表明对任意步长h,梯形法都是绝对稳定的。4.5 线性多步法4.5.1线性多步法的一般公式前面给出了求解初值问题(4.1.1)的单步法，其特点是计算时只用到!的值，此时凡仍，必_1必的值均已算出.如果在计算儿+

23、1时除用話的值外，还用到人亠儿亠丿十的值，这就是多步法.若记忑=起+总,人为步长,九=7(%九)3=0丄屈，贝y线性多步法可表示为Jt-1I儿+1=工务加)+曲工30巴讥T北2=02=-1(4.5.1)其中矽妙为常数，若魔-1+011h,称(4.5.1)为线性k步法计算时用到前面已算出的k个值沧儿小J51.当0j=0时，(4.5.1)为显式方法，当0-1工则称(4.5.1)为隐式多步法隐式方法与梯形方法一样，计算时要用迭代法求畀时I多步法(4.5.1)的局部截断误差定义也与单步法类似.定义5.1设y(x)是初值问题(4.1.1)的精确解，线性多步法(4.5.1)在+1处的局部截断误差定义为ft

24、-1ITn+1=XA+1)一工3-J一应工肉”(7)2=02=-1(4.5.2)若Tb+i=0(k?+1),贝y称线性多步法(4.5.1)是p阶的.如果我们希望得到的多步法是p阶的，则可利用Taylor公式展开，将丫1在处展开到妒粮阶，它可表示为陰1=Uq(a)+U附氏)+U/叭兀)+务护1畀网(耳)+0(护3(4.5.3)注意，(4.5.2)式按Taylor展开可得Jt-1k-1Tn+1=M耳+閒-工口(耳-辺)-曲工-认)2=02=-1+”)+歸导(财一工务肌心)一也”(心3=0占-应工04(心)-如()+才时(耳)+刈2=-1经整理比较系数可得(4.5.4)若线性多步法(4.5.1)为p

25、阶，则可令c厂qqpq+jo于是得局部截断误差y、dg+o旷)(4.5.5)右端第一项称为局部截断误差主项.5+1称为误差常数要使多步法(4.5.1)逼近初值问题(4.1.1),方法的阶pl,当p=1时，则0-G-0,由(4.5.4)得k-1k-1%+呦+-+塔1=1，工加i一工肉=T2-121称为相容性条件.公式(4.5.1)当k=1时即为单步法，若0-1=,由(4.5.6)则得式(4.5.1)就是兀+1=儿+虬，即为Euler法.此时6=0,方法为p=1阶.若得0-1=,由G得I为确定-1及00,必须令1=；=,由(4.5.4)得0一1+炕=1此时(4.5.1)就是儿+1=儿+(+1)即为

26、梯形法.由=4一4久冷一卜一右Z+i一存卩5)+(小故p=2,方法是二阶的,与4.1节中给出的结果相同.实际上，当k给定后，则可利用(4.5.4)求出公式(4.5.1)中的系数口i及以,并求得匚+1的表达式(4.5.5).4.5.2 Adams显式与隐式方法形如Jt-iyn+i=儿+应工0汎斗氏T，匕(4.5.7)的k步法称为Adams方法，当一时为Adams显式方法，当工时，称为Adams隐式方法.对初值问题(4.1.1)的方程两端从耳到耳+1积分得y（耳+i）-yM=o4HK兀E伽显然只要对右端的积分用插值求积公式，求积节点取为G-Z05即可推出形如(4.5.4)的多步法，但这里我们仍采用

27、Taylor展开的方法直接确定(4.5.4)的系数魚e=T,0,用-1).对比(4.5.1)可知，此时%=5=也=%=，只要确定9_屆,，质即可.现在若k=4且0-1=,即为4步的Adams显式方法=1-1-=0儿+1=儿+曲(0。托+01兀一1+炖一2+其中队用a矢为待定参数，若直接用(4.5.4)，可知此时厂一自然成立，再令C=S=5=5=0可得刃+屈+炖+0厂1炖+2炖+加厂冷01+4炖+%冷炖+E：炖+27=-解此方程组得由此得到于是得到四阶Adams显式方法及其余项为(4.5.8)(4.5.9)若0-1工,则可得到p=4的Adams隐式公式，则k=3并令由(4.5.4)可得91951

28、解得乩厂习屁=习屁一帀炖=牙,19720于是得到四阶Adams隐式方法及余项为(4.5.10)(4.5.11)一般情形，k步Adams显式方法是k阶的，k=1即为Euler法，k=2为k=3时，k步隐式方法是(k+1)阶公式，k=1为梯形法，k=2为三阶隐式Adams公式k步的Adams方法计算时必须先用其他方法求出前面k个初值戸仇-)才能按给定公式算出后面各点的值，它每步只需计算一个新的f值，计算量少，但改变步长时前面的儿-1风亠丁2+1也要跟着重算，不如单步法简便.例4.6用四阶显式Adams方法及四阶隐式Adams方法解初值问题y=-+x+!?ox1丿(0)=1,步长h=o.i用到的初始

29、值由精确解y(或f+卞计算得到.解本题直接由公式(4.5.8)及(4.5.10)计算得到.对于显式方法，将=-y+x+1直接代入式(4.5.8)得到对于隐式方法，由式(4.5.10)可得到人+i=片+吕一人+1+耳+i+1)+-5扎_、+血直接求出人+i,而不用迭代，得到几+1=善几+況沁+i+9+19/b-镁_i=2;3?-;9计算结果如表所示.鬲精确解y(跖)=疋呦+鸯.Adams显式方法Adams隐式方法阳|y（跖）涵|%|y（鬲）负10.31.040818221.04081S012.1X1070.41.070320051.070322922.87X1061.070319663.9X10

30、70.51.106530661.106535484.82xl0-fi1.106530145.2X1070.61.148811641.148818416.77X1061.148811016.3X1070.71.193585301.196593408.10X1061.1兀站4597.1X1070.81.249328961.249338169.20xl0-fi1.249328197.7X1070.91.306569661.306572g.g&xiw1.306568848.2X1071.01.367879441.36788P961.05X1061.3678785P8.5X1074.5.3 Adams预

31、测-校正方法上述给出的Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步要做迭代，计算不方便为了避免迭代，通常可将同阶的显式Adams方法与隐式Adams方法结合，组成预测-校正方法以四阶方法为例，可用显式方法(4.5.8)计算初始近似7，这个步骤称为预测(Predictor),以P表示，接着计算f值(Evaluation),这个步骤用E表示，然后用隐式公式(4.5.10)计算+1，称为校正(Corrector),以C表示,最后再计算R+1=了氏小儿+1),为下一步计算做准备整个算法如下：预测：P：灯儿+善(包-59/m_1+37/b_2-Vm_3)求值：E：必=/(兀+if+

32、i)亠亠十血(4.5.12)枚正C：儿+】=儿+可9打+19兀一”门+j求值丑：+1=/(耳+i，片+i)公式(4.5.12)称为四阶Adams预测-校正方法(PECE).利用(4.5.8)和(4.5.10)的局部截断误差(4.5.9)和(4.5.11)可对预测-校正方法(4.5.12)进行修改，在(4.5.12)中的步骤P有对于步骤C有列耳+1)-冗+1盏应夕讥兀)两式相减可得于是有7仪曲)-Yn+1-yLi)若用X+i代替上式-并令显然儿+1,儿+1比唸纭更好，但注意至卩汙的表达式中洛1是未知的，因此改为赠5+】+挽宜-加）F面给出修正的预测-校正格式（PMECME）.+厲-为人+刀九厂汀

33、心）：y=球i+；一心u:冗+】=儿+鲁曲+1汀厂仏+心）(4.5.13)m+iym+i270+1ym+i)应：+i=儿+i)经过修正后的PMECME格式比原来PECE格式提高一阶.4.5.4 Milne方法与Hamming方法与Adams显式方法不同的另一类四阶显式方法的计算公式形如(4.5.14)几+i=尹心+城弘+AZ-i+炖兀-2）这里炕，01，”2为待定常数，此公式也是k=4步方法，即计算儿+1时要用到儿几7儿-沙z4个值为了确定炖個，炖，当然可以利用公式（4.5.4）直接算出，但下面我们直接利用Taylor展开式确定,力。，矗，“2,使它的阶尽量高.方法（4.5.14）的局部截断误

34、差为Tm+1=尹（+血）+巩-陶+却炕卩（耳）+（耳-曲）+02卩（耳-以）.将它在禺点展成Taylor级数，得7JT-7-JTI+1=y（兀）+妙（兀）+药八兀）+页”“（耳）+石汕怯）+側（兀）+-血）-銅g+竽”g-粤才g+粵匕）1譬g+咻心）+加仇）-如乜）+少乜）-分）+普严氏）-.+炖；仏）-2叨g+警严氏）-警尹怯）+畔产牝）-=1+3-（炕+01+02）險（兀）+（|-|+A+202）沪”g+G+务訥-細曲g+G嗚+协+詛約）+学訥岭炖池）+w要使公式的阶尽量高，要令前3项系数为0.即4-凤+01+02)=0厂4+0+幼2=+402=手的系数为0,故848解得A=-亍，0o=三

35、代入公式,(4.5.15)于是得四阶方法(4.5.16)称为Milne公式，它的局部截断误差为（4.5.15）.与（4.5.16）配对的隐式方法为k=3的多步法，它的一般形式可表示为Yn+1=術儿+眄儿-1+勺儿-2+城0-1+1+必兀+01-1）要求公式的阶p=4，可直接用（4.5.4），并令0=6=67=匚=0,可得a0+a1+a2=1_眄_2a2+0_i+A,+0i=1彳a1+4a2+20_-20】=1(4.5.17)-鬥-Sa2+30_+=1鬥+16a2+40_i-40i=1若令也=,可解出筍=0卫1=l,0_i=尿冷炕冷,于是得到下列四阶方法(4.5.18)称为Simpson公式,它

36、的局部截断误差为几+1=-存为归)(耳)+0(胪)(4.5.19)用Simpson公式与Milne公式(4.5.16)相匹配，用(4.5.16)做预测，(4.5.18)做校正，由于(4.5.18)的稳定性较差,因此通常较少使用.为了改善稳定性,可重新选择四阶的隐式公式，Hamming通过试验，发现在(4.5.14)中若令內一,得到的公式稳定性较好，此时(4.5.14)的解为气現宀=-善仏气他兔屁=于是得四阶多步法儿+1=(映一兀-+y(+i+2-Z-1)(4.5.20)称为Hamming公式，它的局部截断误差为几+L-占泊严)+。炉)(4.5.21)用Milne公式(4.5.16)与Hammi

37、ng公式(4.5.20)相匹配，并利用截断误差公式(4.5.15)与(4.5.21)改进计算结果类似Adams预测-校正格式(4.5.13)，可得以下的预测-校正格式(PMECME)：(4.5.22)例4.4用四步四阶显式Milne公式及三步四阶隐式Hamming公式解初值问题y=-+x+i;oxi(o)=i,步长h=o.i初值乃仍由精确解池SX给出，要求计算到冷=二为止，给出计算结果及误差，并与例4.6结果比较.解直接用公式(4.5.16)及(4.5.20)计算.用Milne法计算公式为兔+1=儿+X(須一_1+2人4=3,4其中気=_”+召+1血=_尹+為+=0f1=y1+x1+=-1.0

38、0483742+1.1=0.09516258f2=-y2+x2+1=-1.01873075+1.2=0.18126925/3=-乃+花+1=-1.04081822+1.3=0.25918178九=丹+x(2/3-/2+2/)=1.07032260兀=一兀+可+1=0.32967740旳=戸+罕況(2/4-/3+2/2)=1.10653229误差|y(x4)-y4|=2.55xW6;y(x5)-兀卜1.63x10用Hamming方法(4.5.20)计算公式为可解得n=2,3,4丹=右心叽-兀）+黯*+1+2広-启）=Xxg.168576+xl.56737592=1.040818028.38.3=

39、_乃+也=0.25918198几=存（观-戸）+护（习+1+豪-血=1.皿1颁/4=-儿+心+1=0.32968034k存（妝-E+挣色+1+2了O=1.W653010误差显）一划=2.0xl0_y（x4）-风|=3.9xl0_7y（x5）-5|=5.610-7从所得结果可见Milne方法误差比显式Adams方法误差略小，而Hamming方法与隐式Adams方法误差相当.例4.8将例4.4的初值问题用修正的Milne-Hamming预测-校正公式计算5及兀，初值兀，戸，兀，乃仍用已算出的精确解，即兀=1,”=1.00483742,jp2=1.01873075,73=1.04081822，给出计

40、算结果及误差解根据修正的Milne-Hamming预测-校正公式(4.5.22)得yl=1.106532997=7/+罟心兀一兄）T10小和64了也汁）=y+巧+1=0.393469636X=卜（独一旳）+耳了厲带）+2/4-/3=1.106530419旳=加一善纵応一卅）=1.10653061孔=一旳+花+1=0.39346939y=y,+罟2虫-+2Q=1.14881135膚=-尹譽+阮+1=0.45118865冗=卜（观一幵）+育+2虫/4）=1.14881144兀=冗一善冥舵一出）=1.14閔1161UE）-兀1=4.97x0;|X）-1=2.61x10从结果看，此方法误差比四阶Adams隐式法和四阶Hamming方法小，这与理论分析一致.讲解：线性多步法（4.5.1）的局部截断误差定义为与单步法相似，可表示为（4.5.2），即ft-1A.-1瞪+1=（+却-另-认）-曲另0护（兀-也）2=02=-1只要直接将右端各项在兀处展成Taylor公式，根据公式阶数为西阶，即爲+1=认两）按血的幕整理，令卅卅府各项系数为0则可求得相应的线性多步法及其局部截断误差，这里只用到一元函数的Taylor展开。因此不必记系数满足的公式（4.5.4）,只要直接展开即可，它不但可以求出Adams显式与隐式公式以及Milne公式，H