整理常微分方程数值解

1、精品文档第七章常微分方程数值解7.1 引言本章讨论常微分方程初值问题的数值解法，这也是科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法通常我们假定（7.1.1）中f（x,y）对y满足LipscMtz条件，即存在常数L0,使对内1，gE&有(“1)-14回-R(7.1.2)则初值问题（7.1.1）的解存在唯一.假定（7.1.1）的精确解为尸5）,求它的数值解就是要在区间【上的一组离散点仆=而为由上求尸的近似孙乃,，八.通常取）,h称为步长，求（7.1.1）的数值解是按节点/勺=12）的顺序逐步推进求得巧?.首先，要对方程做离散

2、逼近，求出数值解的公式，再研究公式的局部截断误差，计算稳定性以及数值解的收敛性与整体误差等问题.7.2 简单的单步法及基本概念7.2.1 Euler法、后退Euler法与梯形法求初值问题（7.1.1）的一种最简单方法是将节点二的导数/（/）用差商代替，于是（7.1.1）的方程可近似写成(7.2.1)*y(xn)+吹%,yX4),口=0,L从与出发尸贝湎）一%,由（7.2.1）求得M/）即，Q皆（而3）=乃再将乃指夕（勺）代入（7.2.1）右端，得到”河）的近似乃=修+研勺5）,一般写成%+1=乂+姓（/jj福=。产.(7.2.2)称为解初值问题的Euler法.精品文档精品文档Euler法的几何

3、意义如图7-1所示.初值问题(7.1.1)的解曲线y=y(x)过点P式殉，尸”，从而出发，以产(沟，为)为斜率作一段直线，与直线*=勺交点于耳(看,修),显然有泗=儿十型。，,。)，再从写出发，以八的，力)为斜率作直线推进到上一点斗，其余类推，这样得到解曲线的一条近似曲线，它就是折线此片舄.Euler法也可利用W气+】)的Taylor展开式得到，由-：-)上一工二；十一1:、_：、二。./，)(7.2.3)略去余项，以为、道)，就得到近似计算公式(7.2.2).另外，还可对(7.1.1)的方程两端由二到羽川积分得.|1.：-:,(7.2.4)勺若右端积分用左矩形公式，用心“M&),外+1双双瑞

4、+】)，则得(7.2.2).如果在(7.2.4)的积分中用右矩形公式，则得(7.2.5)精品文档精品文档称为后退(隐式)Euler法.若在(7.2.4)的积分中用梯形公式，则得jiij-,二)I八4丁二一J-U-(7.2.6)称为梯形方法.上述三个公式(7.2.2),(7.2.5)及(7.2.6)都是由几计算八十】，这种只用前一步即可算出外川的公式称为单步法，其中(7.2.2)可由为逐次求出巧小如的值，称为显式方法，而(7.2.5)及(7.2.6)右端含有力匹由八川)当f对y非线性时它不能直接求出八此时应把它看作一个方程，求解八川，这类方法称为稳式方法.此时可将(7.2.5)或(7.2.6)写

5、成不动点形式的方程yn+i=+这里对式(7.2.5)有/L瑞，对(7.2.6)则/=居十号/人，八)，g与八川无关，可构造迭代法匚二：h：(727)由于力X，)对y满足条件(7.1.2)，故有卜储rJ明/a用必1)一/X+Qz)M歌斓-八当自班1或也专，迭代法a?.7)收敛到八一，因此只要步长h足够小，就可保证迭代(7.2.7)收敛.对后退Euler法(7.2.5)，当h；时迭代收敛，对梯形JLr法(7.2.6)，当宓7时迭代序列收敛.例7.1用Euler法、Bt式Euler法、梯形法解精品文档精品文档y|=-y+K+Ly(0)=1取h=0.1,计算到x=0.5,并与精确解比较.解本题可直接用

6、给出公式计算.由于K心力=一丫+工+1$=01礴=0,尤=1,Euler法的计算公式为rs+i=s+K-yx+/+1)=(1-增八+比+血=0-9办+01/+0.1n=0时,外0一9乂+0.H+0./1.n=1,2,3,4的计算结果见表7-1.对隐式Euler法，计算公式为=Lm+A(-A+l十五”1+1)解出八77+砧+1十/=AM十oX十011)1+1.1当n=0时，m=白(为+0为+0D=1P09091.其余n=1,2,3,4的计算结果见表7-1.表7-1例7.1的三种方法及精确解的计算结果段Euler法yn隐式Euler法yn梯形法为精确解贝4)01111011.0000001.00Q

7、0311.0047621.0048370.21.010000L0264461.0135941.019731031.029000L0513151.0406331040818041.05610010330141070097L0703200.5L09D4901.12092211052781.106531对梯形法，计算公式为精品文档精品文档y/i=y”+a+5+i)+(一乂乜+D解得-tt(2一居+A区+&+*)+2A2+=。一次+0一2/+0-21)w11.当n=0时，心一。.9+0.21)004762.其余n=1,2,3,4的计算结果见表2ali7-1.本题的精确解为咫表7-1列出三种方法及精确解

8、的计算结果.7.2.2单步法的局部截断误差解初值问题(7.1.1)的单步法可表示为一(7.2.8)其中与有关，称为增量函数，当含有+】时，是隐式单步法，如(7.2.5)及(7.2.6)均为隐式单步法，而当不含八*1时，则为显式单步法，它表示为L计=+我病(工2L1近)，用=(7.2.9)如Euler法(7.2.2),犬工汽坨=%工力.为讨论方便，我们只对显式单步法(7.2.9)给出局部截断误差概念.定义2.1设y(x)是初值问题(7.1.1)的精确解，记Tn+1“/)一双/)-力贝与，町(7.2.10)称为显式单步法(7.2.9)在/*1的局部截断误差.精品文档精品文档Tz之所以称为局部截断误

9、差，可理解为用公式（7.2.9）计算时，前面各步都没有误差，即丫乳=我时），只考虑由%计算到不用这一步的误差，此时由（7.2.10）有1M步氏，乂-M“+i）一尸（/卜为贝弓乂/）,前=北.局部截断误差（7.2.10）实际上是将精确解式以代入（7.2.9）产生的公式误差，利用Taylor展开式可得到二.e）.例如对Euler法（7.2.2）有奴工,y&）二丁住）,故4必=同工1）-兄4）-好，尸（/）道演+尬-丁区）-犷区J=：力3（小）+白？以匹）十=0（后）它表明Euler法（7.2.2）的局部截断误差为，）+。四,称为局部截断误差主项.定义2.2设式,）是初值问题（7.1.1）的精确解，

10、若显式单步法（7.2.9）的局部截断误差歹是展开式的最大整数，称中为单步法（7.2.9）的阶，含的项称为局部截断误差主项.根据定义，Euler法（7.2.2）中的。=1故此方法为一阶方法.对隐式单步法（7.2.8）也可类似求其局部截断误差和阶，如对后退Euler法5 7.2.5）有局部截断误差精品文档精品文档Zu=J-Gpg)-丘，=yg+&)一双)-砒/+A)=如%)+g-)+-用火/)+刎)+=一?/&)+5/)故此方法的局部截断误差主项为-5墨)亦=1,也是一阶方法.对梯形法26 7.2.6)同样有hT*广以仆+1)7氏)一+/(XqMk)=yK+)-X)-|XUJ+?(+孙=-舐气)+

11、。)它的局部误差主项为与尸Op=2,方法是二阶的.XW7 .2.3改进Euler法上述三种简单的单步法中，梯形法(7.2.6)为二阶方法，且局部截断误差最小，但方法是隐式的，计算要用迭代法.为避免迭代，可先用Euler法计算出外+i的近似无+】，将(7.2.6)改为艮+1=八+舒(/必)，入_九川二居十万【/1赤，居)+/(演+卜见七1)1(7.2.11)称为改进Euler法，它实际上是显式方法.即h1-1一二(7.2.12)吕右端已不含%-1.可以证明T*+=Og,声=2,故方法仍为二阶的，与梯形法一样，但用(7.2.11)计算%+1不用迭代.精品文档精品文档例7.2用改进Euler法求例7

12、.1的初值问题并与Euler法和梯形法比较误差的大小.解将改进Euler法用于例7.1的计算公式1=%+1d*+1）+&+D+x+biria（2-A）ik（y-h）hfj肛27）、=1-氏+5（/+力+?）乙乙qda-=0905八+00955+01当n=0时，%=900+。95%+01-1一005000.其余结果见表7-2.表7-2改进Euler法及三种方法的误差比较A改进E必升法打误差|外.烬乱）|Eu加方法-婷）|梯形法1%-加力0.110050001.6X10448xl0-375x100.2L01902529x1。-*87xl0-31.4x10031.0412184.0xl0M1.2x1

13、0“19x1。-*0.410708024.8x10m1.4xlO-32.2X10-411070765,5x10-41.6xl0-a2.5*10-从表7-2中看到改进Euler法的误差数量级与梯形法大致相同，而比Euler法小得多，它优于Euler法.讲解:求初值问题（7.1.1）的数值解就是在假定初值问题解存在唯一的前提下在给定区间加上的一组离散点值;而网父父/）-加（弓“/）=。（小】），定义局部截断误差，它表示用精品文档精品文档精确解代入计算公式（7.2.9）产生的公式误差为（圾“）田越大表明公式逼近微分方程的精度越高，因此就定义为公式的阶，通常尸之1的公式才能用于计算初值问题（7.1.1

14、）的数值解。利用Taylor展开时，只要将看+1的表达式在右处展开成Taylor公式就可得到不同公式的局部截断误差。如7.2.2所给出的Euler法。后退Euler法和梯形法，它们只需用一元函数的Taylor展开，与后面7.5节的多步法完全一致，而通常单步法（7.2.9）的一般情况则需要用二元函数的Taylor展开，才能得到公式的具体形式和局部截断误差。例如对改进Euler法,其局部截断误差由（7.2.12）可得%=JM4）-y/氏，以勺）+砂（仆”勺）要求出它的结果就要用到二元函数兀力的Taylor展开，将在7.3节再作介绍。7.3Runge-Kutta方法7.3.1 显式Runge-Kut

15、ta法的一般形式上节已给出与初值问题（7.1.1）等价的积分形式(7.3.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题（7.1.1）的数值方法，若用显式单步法%+】=%十娜（冬以E=。工(7.3.2)当=胴.,%），即数值求积用左矩形公式，它就是Euler法（7.2.2）,方法只有一阶，若取(7.3.3)步区,以）=5芍,八）+/0用入十财区,八）精品文档精品文档就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(7.3.1)右端积分表示为f(

16、x,y(x)dx=筏&及yg+白河)+“(犷*】)扁i-i注意，右端f中产&n+%h)还不能直接得到，需要像改进Euler法(7.2.11)一样，用前面已算得的f值表示为(7.3.3)，一般情况可将(7.3.2)的表示为r，二二，；(7.3.4)其中;11川口,区，/就A+人工”号)02,3,/)这里c”出晶b-12,rj=12J-D均为待定常数，公式(7.3.2),(7.3.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即，上工)，而ki由前面(i-1)个已算出的孰,曷，,号_表示，故公式是显式的.例如当r=2时，公式可表示为一：(7.3.5)其中用八)+%

17、八十瓯岫).改进Euler法(7.2.11)就是一个二级显式R-K方法.参数“.右町&】取不同的值，可得到不同公式.7.3.2二、三级显式R-K方法对r=2的显式R-K方法(7.3.5),要求选择参数也使公式的阶p尽量高，由局部截断误差定义.一；、.+1(7.3.6)精品文档精品文档令&=双小）,对（7.3.6）式在（私，人）处按Taylor公式展开，由于h2h34y（/G=yg+h）=y”十如十方国尸十伙期）Z=/（/，片）=力d=瓦/区t八）=/X，八）+g，居，（小+出居+%竭户/氏,八）+&又）m臬）（域i）+。5）将上述结果代入（7.3.6）得*,.X也十万&K）十力氏，八）十*）-

18、Me遥+的a+M区.十匹（仆区）-。）力=（1-6-q）妨十gr勺勺百/（4,）+一心次）为区以）+0（*）要使公式（7.3.5）具有的阶p=2,即T/10（/），必须11-二-.-（7.3.7）口L.T1,1即：一由此三式求犯当】的解不唯一.因r=2,故工/口,于是有解1-一一,-丁（7.3.8）它表明使（7.3.5）具有二阶的方法很多，只要勺k都可得到二阶R-K方法.若取:S=：,则与=；，%=%i=1,则得改进Euler法（7.2.11）,若取9=1,Zj乙则得=。色=%1=；,止匕时（7.3.5）为W-精品文档精品文档:一.；.、.(7.3.9)其中:称为中点公式.后退Euler法(7

19、.2.11)及中点公式(7.3.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(7.3.6)的展开式中要增加3项，即增加三个方程.加上(7.3.7)的三个方程求4个待定参数是无解的.当然r=2,p=2的R-K方法(7.3.5)当0Ho取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即照一y2h)=c内+%上”6&其中:i二Y二%三/g+%瓦y*+与/也+%麻1将七段按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得1+1=求/+旬以妣/，外出=5犷)可得三阶

20、方法，其系数应满足方程2%=6这是8个未知数6个方程的方程组，解也是不唯一的，通常右0.一种常见的三级三阶R-K方法是下面的Kutta三阶方法:精品文档精品文档一二-二（7.3.11）口用=/出）描三/X+/，乂+：用）用=（4+瓦%-%+2力尼）7.3.3四阶R-K方法及步长的自动选择利用二元函数Taylor展开式可以确定（7.3.4）中r=4,p=4的R-K方法，经典的四阶R-K方法是：二十一二十幻一二二一二%+，：（7.3.12）6月=0%十；及/性）网=/（/+卜+;3）一（“十机M+蚂）它的局部截断误差%二承、故p=4,这是最常用的四阶R-K方法，数学库中都有用此方法求解初值问题的软

21、件.这种方法的优点是精度较高，缺点是每步要算4个右端函数值，计算量较大.例7.3用经典四阶R-K方法解例7.1的初值问题尸十仍取h=0.1,计算到株=0-5,并与改进Euler法、梯形法在知=。-5处比较其误差大小.解用四阶R-K方法公式（7.3.12）,此处而=0,%=1八01于是当口=。时精品文档精品文档孔-/（即丹）=/+1=0/7（配+；及必+、D乙M=（-1+1龙+。-3）。+1=005电二丁（两+g机4+玲）二（7+一.+3九+（喂+学两+10375&=/（%+机汽+为）二M冷h1”=（_1+）为+（17+）/+1=009525乙乙I,于是冯=为+打电+2尾+2&+质）=1+xO.

22、29025=1.00483750,按公式66（7.3.12）可算出乃=1018730903=1.040818422=1,07032029,=11065309此方法误差：-1.:1改进Euler法误差：1-1：梯形法误差：,1：-可见四阶R-K方法的精度比二阶方法高得多.用四阶R-K方法求解初值问题（7.1.1）精度较高，但要从理论上给出误差限哦嘘的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h?通常是通过不同步长在计算机上的计算结果近似估计.设双仆）在飞处的值外.尸（/），当匹+冬会时，#湎）的近似为喏i,于是由四阶R-K方法有精品文档精品文档若以:为步长，计算两步到则

23、有JL次（2）于是得ym+i.?）Mx）一、&即：】.:二一.【二：或陶1（-）R武/+1）-强*lAii-XUI（7.3.13）它给出了误差的近似估计.如果乙，匕：-y驾I，义为给定精度），则认为Jfc,/j以5为步长的计算结果满足精度要求，若记1焉一斓1宾,则还可放大步长.因此（7.3.13）提供了自动选择步长的方法讲解：求初值问题（7.1.1）的单步法主要是指Runge-Kutta法，本节主要讨论显式RK方法，建立具体的计算公式使用的是Taylor展开，形如（7.3.4）的显式RK方法，当r=1时就是Euler法，因此只要讨论,，之2的计算公式，在r确定后如何推导公式都是一样的，只是r越

24、大计算越复杂，为了掌握了解公式来源，只要以r=2为例推导计算公式即可。因此本节重点就是用Taylor展开求出r=2的显式R-K方法的计算公式，由于方法的局部截断误差为（7.3.6）,4+1的右端有/+%氏+%】）的项，要对它做Taylor展开，就要用到二元函数克J）的Taylor展开，按照二元函数Taylor级数精品文档精品文档/(x+Ax,y+Ay)=+0(hr+1)2团改办JV、本（了J）,亨（K，T）11r/Jl心不居尸）=/）+的-+7713工+(7.3.14)由吻2!微2&xLy3兀疗）+（4了觉八支）+将它用到（7.3.6）的+i的展开式中，即可得到按人开幕整理出的结果，对r=2的

25、公式只能得到P=2阶的公式，即工+L。仍与，于是2级R-K方法（7.3.5）的系数，与，的法?1必须满足（7.3.7）给出的方程，它的解由（7.3.8）给出，只要门M求出的公式都是r=2的2阶R-K方法。而常用的就是二得到的2改进Euler法（7.2.11）和口二1得到的中点公式（7.3.9）。7.4单步法的收敛性与绝对稳定性7.4.1 单步法的收敛性定义4.1设y（x）是初值问题（7.1.1）的精确解，乃是单步法（7.3.2）在%=%+泌处产生的近似解，若则称方法（7.3.2）产生的数值解收敛于WQ.实际上，定义中，是一固定点，当h一0时n-8,n不是固定的.因由=福显然方法收敛，则在固定点

26、工.处的整体误差/=尸（/）-八=03,），当pi1M_A时TF.下面定理给出方法（7.3.2）收敛的条件.定理4.1设初值问题（7.1.1）的单步法（7.3.2）是p阶方法（p1）,且函数对y满足Lipschitz条件，即存在常数L0,使对的衣仁衣，均有|灰%y（h）-,h）|I向1=A,显然计算是不稳定的.如果用后退Euler法(7.2.5)解此例，仍取h=0.025,则儿期几十灰+】=乂-2,5八+】，即八.1=与乂Him显然当,0,计算是稳定的.由此看到稳定性与方法有关，也与力有关，在此例中力=-10.在研究方法的稳定性时，通常不必对一般的f(x,y)进行讨论，而只针对模型方程精品文档

27、精品文档(7.4.2)y=fRZ)0这里尤可能为复数.规定Re&)是因为RK。时微分方程(7.4.2)本身是不稳定的，而讨论数值方法(7.3.2)的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题(7.1.1),若将f(x,y)在处线性展开，可得出尤y)七f(/，然)+；(公,八)&-/)+4区必)(y-八)于是方程(7.1.1)可近似表示为9=由+E&型乱)，兄=f；炽n，%)它表明用模型方程(7.4.2)是合理的，至于模型方程(7.4.2)中所以用复数入是因为初值问题(7.1.1)如果是方程组，即蚱E，曰立二则是(mxm)阶矩阵，其特征值可能是复数.当然对单个方程，入就是实

28、数，此时只要规定“+1（M）3+A（m）定义4.2将单步法（7.3.2）用于解模型方程（7.4.2），若得到（7.4.3）中的|E（h|1则称方法是绝对稳定的.在复平面上复变量满足E&兄）|41的区域，称为方法（7.3.2）的绝对稳定域，它与实轴的交点称为绝对稳定区间.例如对Euler法，1磔兄）1=|1+卜兄|41在复平面”,上是以（-1,0）为圆心,以1为半径的单位圆域内部，当片为实数时，则得绝对稳定区间为-220,因A0,故有0h在例7.4中4心_=0.02时方法稳定，而例中取.；h=0.025故不稳定.对后退Euler法（7.2.5）,匕+】=居+%+1因片0,故卜必1l”，其绝对稳定

29、域是以（1,0）为圆心的单位圆外部，绝对稳定区间为-8（应之0方法都是绝对稳定的.二阶R-K方法的绝对稳定区间为-2A20.三阶R-K方法的绝对稳定区间为-2一51税0.四阶R-K方法的绝对稳定区间为-2785根0.例7.5用经典四阶R-K方法计算初值问题y=-2oxoi）（o）=1步长取h=0.1及0.2，给出计算误差并分析其稳定性.精品文档精品文档解本题直接按R-K方法（7.3.12）的公式计算.因精确解为/=，其计算误差I%-如表所示.020406031.0h-0.10.0927950.012010000136C00001520.000017h=0.24g825.C125.0625.03

30、125.0从计算结果看到，h=0.2时误差很大，这是由于在入=-20,h=0.2时入h=-4,而四阶R-K方法的绝对稳定区间为-2.785,0，故h=0.2时计算不稳定，误差很大.而h=0.1时h；l=-2,其值在绝对稳定区间-2.785,0内，计算稳定，故结果是可靠的.讲解：由于微分方程初值问题数值解公式求出的解酒外，八是一个逐次递推的过程，因此原始数据误差及计算过程舍入误差对解的影响就是数值方法绝对稳定性研究的问题，如果由人计算尸z误差不增长，方法就是绝对稳定的。为使问题得到简化通常就是将方法用于解模型方程（7.4.2）,对于单步法得到的差分方程为二田:山K,由于模型方程的冗川二廿,代入E

31、uler法+L八+灰=乂十研右二。十娘）久，得物）=1+麻，对二阶R-K方法，例如，用改进Euler法%+1=+维丁（羽，居）+/（%+比其+&，八于是=+-+Aa+)=豆（瓦I）=1+冉几+1侬冗）2101M伽。=1+&丸牛十上（蚁广对三阶R-K方法有23!精品文档精品文档5(2)=1+2+13建9+-(hx?+(烟对四阶R-K方法有23!41只要冏陋义)卜1方法，就是绝对稳定的，这时上的值当n增大式是减少的,故计算稳定。这时舍入误差影响可忽略不计，而当忸5上1,则卜增大，方法不稳定，计算结果是不可靠的。因此用显式单步法必须使阿砌卜1,也就是步长选择要满足这一要求。对于隐式的梯形公式一【+：

32、h将模型方程上力，即二妙代入得瑞十小八十11+工义兀_2于是八十1=-i-八注意4(/)+-必卜(/)t&k)+ftyaJ+1i=0dffi血)不-人工43)-Q(/)+i=-l经整理比较系数可得精品文档精品文档Jt-1i=f)i-1Jt-1，G=一(7)%-工兑(7.5.4)S=0i=-lIE-1。广产学*-汨自F”尸g,白若线性多步法(7.5.1)为p阶，则可令1=,=C,=O，CiHO于是得局部截断误差二，一：一一(7.5.5)右端第一项称为局部截断误差主项.Cmi称为误差常数.要使多步法(7.5.1)逼近初值问题(7.1.1),方法的阶p1,当p=1时，则Co=g=O,由(7.5.4)

33、得juiT斯+苗4-L%一用=72V111称为相容性条件.公式(7.5.1)当k=1时即为单步法，若尸-1=,由(7.5.6)则得式(7.5.1)就是了.”=+权”,即为Euler法.止匕时Cn=-0,方法为p=12阶.若处=0,由Cl得=1,为确定外及为，必须令G=q=o,由(7.5.4)得a十瓯=1及/5精品文档精品文档此时（7.5.1）就是尸甘一=K十女工+九+1）即为梯形法.,由一一，:-一一_31上！041212故p=2,方法是二阶的，与7.1节中给出的结果相同.实际上，当k给定后，则可利用（7.5.4）求出公式（7.5.1）中的系数的及以，并求得T-i的表达式（7.5.5）.7.5

34、.2Adams显式与隐式方法形如口】(7.5.7)=乂+力7，制=止1戊一,i=l的k步法称为Adams方法，当。时为Adams显式方法，当员】工。时,称为Adams隐式方法.对初值问题（7.1.1）的方程两端从为到积分得y（Kn+l）ryOn）=小巧（工）占了显然只要对右端的积分用插值求积公式，求积节点取为品-叼仆工2即可推出形如（7.5.7）的多步法，但这里我们仍采用Taylor展开的方法直接确定（7.5.7）的系数.G=一2,7）.对比（7.5.1）可知，此时4-0,只要确定员1，尸如，尸X即可.现在若k=4且a1=,即为4步的Adams显式方法儿+1=怎+双凤九+回兀.1十也广一十四工

35、与）精品文档精品文档其中乩,尸卜为，区为待定参数，若直接用（7.5.4）,可知此时ci-%=1-1-o自然成立，再令c1=cq=c4=o可得凤-22-3伪*1人+&M+9区-凤-8凡+27仇解此方程组得由此得到55379乩二函）五屎一正.广彳一不）小行于是得到四阶Adams显式方法及其余项为必+广九+白-59工373TxQ(7.5.8)。卷1g+加）(7.5.9)若瓦1”。，则可得到p=4的AdamS急式公式，则k=3并令G=CCC”0由（7.5.4）可得_1+1+/+同=1B-B三十月十%二；-A-8A=1精品文档精品文档91951113A19解得用西网加血=-M尸不而丁口m”一而于是得到四

36、阶Adams隐式方法及余项为.1.二二二-二一（7.5.10）%一厂二（7.5.11）一般情形，k步Adams显式方法是k阶的，k=1即为Euler法，k=2为kyN+i=A+5。工一工7）Uk=3时，=K+（23力1&力.1+5冗修）.k步隐式方法是（k+1）阶公式，k=1为梯形法，k=2为三阶隐式Adam公式h%,i=%+R5九i+%-九）Xk步的AdamspJ法计算时必须先用其他方法求出前面k个初值u/卜J】）才能按给定公式算出后面各点的值，它每步只需计算一个新的f值，计算量少，但改变步长时前面的片-1，入4】也要跟着重算，不如单步法简便.例7.6用四阶显式Adams方法及四阶隐式Ada

37、ms方法解初值问题产p+工+1,04工,”0）1步长卜=0.1用到的初始值由精确解y=e计算得到.解本题直接由公式（7.5.8）及（7.5.10）计算得到.对于显式方法，将人孤力=-y+耳+1直接代入式（7.5.8）得到%=乂十五6工-59%-9小=34尸乙r精品文档精品文档其中1一-If.对于隐式方法，由式(7.5.10)可得到为+八+学*9居+1+/+i+1)+19%7工.1十九“直接求出八*1,而不用迭代，得到Q1%，=方十浙9%+1+9+19力-5A-1=2,3,，90-324y计算结果如表所示.辱精确解Adams显式方法Adams隐式方法为6K痂凫1为0.31.04081822104

38、0813012.1x10-7041070320051.070322322.S7X1061070319663,9x1W051106530(561106535484.82x10S1.106530N52x10-70.6114S311H1l0-1.148311016.3x10-70.71.193585301.1593408.10x10-511965E4597,1X107081249328961.24的招1(5S.20xl(H1245328197.7x1W09130656966L306579629.P6xlO13065E鸵4SJXIO-T1.01367S7P441.367389

39、96L05X10-5136787059S.5X1Q-77.5.3Adams预测-校正方法上述给出的Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步要做迭代，计算不方便.为了避免迭代，通常可将同阶的显式Adams方法与隐式Adamsf法结合，组成预测-校正方法.以四阶方法为例，可用显式方法(7.5.8)计算初始近似乂雪，这个步骤称为预测(Predictor),以P表示，接着计算f值(Evaluation),丁*三,这个步骤用E表示，然后用隐式公式(7.5.10)计算八+1,称为校正(Corrector),以C表示，最后再计算为下一步计算做准备.整个算法如下：精品文档精品文档h预

40、测：P：匕】=%+五3工-594-1+37c求值*E./4i=/（/+1,*】）1-h（7512）校正C：人+R9公+19工-5九十九）求值距加二（5居+J公式（7.5.12）称为四阶Adams预测-校正方法（PECE）.利用（7.5.8）和（7.5.10）的局部截断误差（7.5.9）和（7.5.11）可对预测-校正方法（7.5.12）进行修改，在（7.5.12）中的步骤P有MlMk）-匕1如丽“叫对于步骤C有1Q-W+-加/y（不）720两式相减可得二二二，L，一于是有251次（%+1-%+】）*-痂居+J193-人】刃-八+D若用瑶口代替上式八A并令濡-丁柒+牛8】丁丁京）居+i-yM+i

41、v(yM+iyM+i)占JW显然尸方以+1比4更好，但注意到尸点的表达式中质1是未知的，因此改为精品文档精品文档了*=泣1+元助;-0)ury卜面给出修正的预测-校正格式（PMECME）.p-yLi=+白(55力-59工37九-9九)，!一1/渭=41+券。:-M)上:篇1n小渭)Cy：+u*s您+19/-5九十九)2419-ym+lJKitzt超一八小1)乙FU产:力+L=/(Xjt+lJjt+1)(7.5.13)经过修正后的PMECME式比原来PEC弗式提高一阶.7.5.4Milne方法与Hamming方法与Adams显式方法不同的另一类四阶显式方法的计算公式形如八川=八与十双凤工+AZ,

42、-1+hH(7.5.14)这里/o,44口为待定常数，此公式也是k=4步方法，即计算八+1时要用到JQVl居也以a4个值.为了确定力。，4,，当然可以利用公式（7.5.4）直接算出,但下面我们直接利用Taylor展开式确定，/。，用伉，使它的阶尽量高.方法（7.5.14）的局部截断误差为T”i十血）十犷氏一生）+从综氏）十&一与十凡尸（乙一物将它在a点展成Taylor级数，得精品文档精品文档%=(&)+A纱O*.)+步区)+-)-湖气)+萼y“(x*)-萼V(%)+粤/)&)-粤*(“)+,】-(h2h31%*A3(/)+屈1/(/)-Q(4)+p/K)-石尸阴(公)+下;(/)-+色LI、n?iizy(2心)(2以)m(2血)*I1产吗)-23口苛*(/)+彳-*(&)-卜L=1+3(其十4+)+gg+同+2)y匹)+12714口、下川/、,18110SaY,4(4j/、(7十下一大片一可儿物尸(/)+(万一否+工用1+/4讷尸(勺)十13s1寸二/一不自3严(QM要使公式的阶尽量高，要令前3项系数为0.即4一