chap7常微分方程数值解ppt课件

1、第七章 常微分方程数值解法1 1 、引言、引言00:( ,)()dyfx yaxbdxyyx一 阶 常 微 分 方 程 的 初 值 问 题2yxy2y4xy4x2yf x y4y 13xdy2yf x y4dxxy 13: - ( , ) ( )-( , )( ) 例例方方程程令令：且且给给出出初初值值就就得得到到一一阶阶常常微微分分方方程程的的初初值值问问题题： ( , )( , )( , )f x yyLipschitzLf x yf x yL yy只要函数适当光滑连续，且关于 满足条件，即存在常数 ，使得由常微分方程理论知，初值问题的解必存在且唯一。 微分方程的数值解：设方程问题的解微分

2、方程的数值解：设方程问题的解y(x)y(x)的存在区的存在区间是间是a,ba,b，令，令a= x0 x1 xn =b,a= x0 x1 xn =b,其中其中hk=xk+1-xk , hk=xk+1-xk , 如是等距节点如是等距节点h=(b-a)/n , hh=(b-a)/n , h称为步长。称为步长。y(x)y(x)的解析表达式不容易得到或根本无法得到，我们用的解析表达式不容易得到或根本无法得到，我们用数值方法求得数值方法求得y(x)y(x)在每个节点在每个节点xkxk上上y(xk)y(xk)的近似值，用的近似值，用ykyk表示，即表示，即yky(xk)yky(xk)，这样，这样y0 , y

3、1 ,.,yny0 , y1 ,.,yn称为微分称为微分方程的数值解。方程的数值解。 主要问题v如何将微分方程离散化，并建立求其数值解的递推公式；v递推公式的局部截断误差，数值解与精确解的误差估计；v递推公式的稳定性与收敛性。v用差商代替微商用差商代替微商v数值积分数值积分vTaylor展开展开微分方程离散化常用方法微分方程离散化常用方法 1,1A ( , ( )nnnnnnnny xy xdyf x y xydx xxx用差商代替微商11111 , , , 0, 1, 2,nnnnnnnnnnnnnnhyyfhxxyxyxyyx yyyxhfny用代替，则 11B. : ( , ) (0,1

4、,) nnnnxxxxdydxf x y dxndx用数值积分方法离散化1111,(), (),( , )(,) (0,1,(,)nnnnnnxnnxnnnnyyy xy xf x y dxhf x yyh xnyfy用代替对右端积分采用取左端点的矩形公式 则有11111111,(), (),( , )(,) (0, 1, )(,)nnnnnnxnnxnnnnyyy xy xf x y dxhf xyyyhf xyn用代替对右端积分采用取右端点的矩形公式 则有11111111,(), (),( , )(,)(,)2 (,)(,)2 (0, 1 ,)nnnnnnxnnnnxnnnnnnyyy x

5、y xhf x y dxf x yf xyhyyf x yf xyn用代替对右端积分采用梯形公式 则有22C ( ) ()()()()2 ()(, ()()2nnnnnnnnnxy xTaylorhy xhy xhy xy xhy xhf xy xy x在附近的展开：11 ( ) () ( ,) 0, 1, 2, nnnnnnnhyy xy xyyhf x ynTaylor取的线性部分， 且得的近似值：展开法不仅可得到求数值解的公式，且容易估计截断误差。1 解常微分方程初值问题的Euler方法EulerEuler方法方法EulerEuler方法的误差分析方法的误差分析v向前向前EulerEul

6、er公式公式EulerEuler折线法或显格式）折线法或显格式）v向后向后EulerEuler公式后退公式后退EulerEuler公式）公式）v梯形公式改进的梯形公式改进的EulerEuler公式）公式）vEulerEuler预估校正格式预估校正格式一、一、EulerEuler方法方法0010, 0 1,1, nnnnny xyyyhfxyxxnhban,NhN1 1、向前、向前EulerEuler公式公式 000000,0000 , , , , x y yy xdyf x yydx xf x yx y几何意义由出发取曲线的切线(存在!)，则斜率由于及已知，必有切线方程。00000000, (

7、) (,)dyyyxxyxxf x yydx x由点斜式写出切线方程：1011000 ,hxxhyyyhf x y等步长为 ，则，可由切线算出：（）11 , 0 1 2 nnnnnyy xxyyhf xyn逐步计算出（ ）在，点值 ：（）,， ，用分段的折线用分段的折线逼近逼近函数逼近逼近函数2、向后后退的、向后后退的Euler 方法方法11 nnny xy xy xh用向后差商：11100, nnnnyyhfxyy xy (隐式算法) 01111E , 0,1,2,1,nnnnkknnnulerhfhfknyyyxyyyx为避免解非线性方程，与法结合迭代法3、梯形公式11 ,nnnnxy x

8、y xfx y dtx由积分途径 11 ,nnnnyy xyy x积分用梯形公式，令则得11100,2nnnnnnhyyfxyfxyy xy(0)1(1)( )111012,2 012nnnnkknnnnnnEulernyyhf x yhyyf x yf xyk, , ,同样与法结合，形成迭代算法，对， ， ，4、改进的尤拉公式、改进的尤拉公式梯形公式虽然提高了精度，但使算法复杂。而在实际计算中只迭代一次，这样建立的预测校正系统称作改进的尤拉公式。1111 ( ,); ( ,)(,),2nnnnnnnnnnyyhf x yhyyf x yf xy预测校正11(,);(,);() / 2.pnn

9、ncnnpnpcyyhfxyyyhfxyyyy二、二、Euler方法的误差分析方法的误差分析11111) () nnnnnnnTy xyyyy xEulery局部截断误差在一步中产生的误差而非累积误差：其中是当（精确解！）时由法求出的值，即无误差！ 1121 () ()()(, () 2nnnnnnnnny xxTaylory xy xhy xhf xy xhyxx将在点展开： 1121111 ,(), 2nnnnnnnnnnnnnyyhf x yyy xhf x y xhTy xyyxx则22212 max( ) , ( ) 2a x bnMy xy xhTMO h 令充分光滑，则： 11

10、, , 1(1)nnnnnnnnnnnnnnnnnyyy xhf x y xyhf x yy xyh f x y xf x yLipschitzhL y xyhL e由条件11000121211 1 N 0 111 1nnnNNNNNNNeThL eney xyeThL eThL ThLThLT对一切 成立,对取定，由， 则：2） 总体方法误差总体方法误差 , ,nnnnnnf xyLipschitzf xy xf xyL y xy递推方法：从任意两相邻步的总体误差关系推出总体误差与步长的关系。由微分方程解的存在唯一性自然假定（ ， ）充分光滑，或满足条件： 1111111111111 1 ,

11、 nnnnnnnnnnnnnnnnnnney xyney xyy xyyyTyyyyyy xhf x第 步的总体截断误差记为则对步 ：以下估计其中 ,ny x221211 11 1NNNNNNTO heThL ThLThLT由局部截断误差，则11200 1 1NNNKkkkkOhLhLhT 211 11NhL O h O hhL000 lim1lim1 NNxxNhhhLhLhLxxeh与 步 长无 关 常 数总体截断误差与局部截断误差的关系是：1O h总体截断误差局部截断误差一般地，方法的总体截断误差阶越高，精度也越高。pO hp定义：一个方法的总体截断误差若为， 则称之为 阶方法。误差分析

12、表误差分析表EulerEuler方法方法局部截局部截断误差断误差总体截总体截断误差断误差迭代收敛迭代收敛条件条件向前向前EulerEuler方法方法O(h2)O(h)向后向后EulerEuler方法方法O(h2)O(h)0hL1梯形公式梯形公式O(h3)O(h2)0hL2(L为为Lip常数）常数）向后向后Euler Euler 方法收敛条件与截方法收敛条件与截断误差断误差 11111111111 , 01kkkknnnnnnkknnh ffhLhLyyyyxxyy收敛条件 21 01 nOO hhLhT局部截断误差，整体截断误差(当时)梯形公式的收敛性梯形公式的收敛性 11111111111

13、,2 2 012kkkknnnnnnkknnhffhLhLyyyyxxyy收敛条件 Euler梯形公式比法的局部与总体误差均高一阶，但每次迭代均多算一次函数值提高精度的计算代价。2(01);(0)1.xyyxyy例：用尤拉公式和改进的尤拉公式解初值问题10.12().nnnnnhxEuleryyh yy解：取步长，公式为： 112();2();1().2npnnnncnppnpcxyyh yyxyyh yyyyy改进的尤拉公式为：20.2(00.6); (0)1.hyyxyxy 例1：取步长，用欧拉法解初值问题122 , 0.2 0.80.2nnnnnnnnnnnEuleryyhfxyyyx

14、yyx y解：格式为：0122 1 0.2 0.8, 0.4 0.6144 0.6 0.461321yyyyyyy由计 算 得0.283(12); (1)2.5hyyxy例2：取步长，用梯形解初值问题小数点后至少保留 位。11111 ,20.2 83832nnnnnnnnnnhyyfxyfxyyyyy解：梯形公式为： 1012345716 131312, 1.22.307691.42.473371.62.562581.82.610622.02.63649nnyyyyyyyyyyyyyy故 由计 算 得 0.23(12); (1)2.5hyxyxy例2 ：取步长，用梯形解初值问题小数点后至少保留

15、 位。111111 ,2 0.1 33nnnnnnnnnnnnhyyfxyfxyyyx yxy解：梯形公式为： 0111110011234511111110.30.1 332,2.6*,&,%,?nnnnkknnnnnnkyyx yyyx yxyyyyyyyyyy迭代格式： 021234522222,*,&,%,yyyyyy20.2sin0 (1)1.1.21.45hyyyxyyy例3：取步长，用欧拉预校方法解初值问题计算及的近似值，小数点后至少保留 位。111121212111,2 0.2sin0.1sin + sinnnnnnnnnnnnnnnnnnnnnnnnyyhfxy

16、hyyfxyfxyyyyyxyyyyxyyx解：欧拉预校格式为： 0112211,0.63171 1.20.7154880.476961.40.52611yyyyyyyy故 由计 算 得 50 (0)1.0.10.210yyyhy例4：写出用反复迭代的欧拉预校法解初值问题的计算公式，并取步长，计算，要求迭代误差不超过。 01111101111,2 0,1,2, 0,1,2, 0.90.950.05nnnnkknnnnnnnnkknnnyyhfxyhyyfxyfxyknfx yyyyyyy 解：欧拉预校的迭代格式为：取 0012111341143751141101, 0.9,0.905,0.90

17、475, 0.9047625,0.9047618756.251010, 0.10.904761875yyyyyyyyyyyy故 由计 算 得由于 是 取 0012232452254652252201, 0.814286,0.818809, 0.818583, 0.818595,0.8185941010, 0.20.818594yyyyyyyyyyyy故 由计 算 得由于 是 取2.2.龙格龙格库塔方法库塔方法 基本思想基本思想 二阶二阶R-K方法方法 三阶三阶R-K方法方法 四阶四阶R-K方法方法 变步长变步长R-K方法方法1112() ()()()()()()2!( )( , ), ( )(

18、 , )( , ) ( , ),nnnppnnnnxypTaylory xyy xhhy xhy xyxyxPy xf x yyxfx yfx y f x y若用 阶多项式近似函数有：其中。但由于公式中各阶偏导数计算复杂，不实用。一、基本思想一、基本思想(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1); ;2,3,jjjjyffffyffxyffyffxyffyffjxy一般地有111112121(,)11()22 (,)(,)nnnnnnnnnnEulerEulerEuleryyhKKf xyEuleryyhKKKf xyKf xh yhK如果将公式与改进公式写成下列形式：

19、公式改进公式11 ( , )()( , )( , )nnf x yy xyf x yf x y以上两组公式都使用函数在某些点上的值的线性组合来计算的近似值。Euler公式：每步计算一次的值，为一阶方法。改进Euler公式：需计算两次的值，二阶方法。( , )(,)( )-nnnf x yxyTaylory xxTaylorR K于是可考虑用函数在若干点上的函数值的线性组合来构造近似公式，构造是要求近似公式在处的展开式与解在 处的展开式的前面几项重合，从而使近似公式达到所需要的阶数。即避免求偏导，又提高了方法的精度，此为方法的基本思想。11111-(,)(,) (2,3,),pnniiinnii

20、ninijjjiijiRKyyhc KKf xyKf xa h yhb Kipabc一般地，方法设近似公式为其中 ，都是参数.,(,)( )iijinnnab cxyTaylory xxTaylor确定参数 ，的原则是:使近似公式在处的展开式与在 处的展开式的前面项尽可能多地重合。二、二阶龙格库塔方法二、二阶龙格库塔方法111221222112()(,)(,)nnnnnnpyyh c Kc KKf xyKf xa h yhb K当时，近似公式为 112221123221( ,)( ,)(, ( ,)( ,) ( ,) ( ,)( ,) ( ,)( )nnnnnnnnnnnnnnnxnnynnn

21、nnx yTayloryyhc f x yc f xa h yhb f x yyh c f x yc f x ya hf x yhb f x y f x yO hy上式在处的展开式为12222321() ( ,)( ,) + ( ,) ( ,)( )nnxnnynnnncc f x y h c a f x yb f x y f x yhO h123123()()()()()()2(,)(,)2 (,) (,)()nnnnnnnnnxnnynnnny xxTaylorhy xy xhy xy xO hhyf xy hfxyfxyf xyO h在 处的展开式为122 232 211 1/2 1/2

22、 ( ),cccacbOh有 无 穷 多 组 解 ， 每 一 组 解 得 一近 似 公 式 ， 局 部 截 断 误 差 均 为这 些 方 法 统 称 二 阶 方 法 。43()()O hO h可以证明，无论这四个参数如何选择，都不能使局部截断误差达到，也即在计算两次函数值的情况下，局部截断误差的阶最高为。122211121211,1,2()/2(,)(,)nnnnnnccabEuleryyh KKKf xyKf xh yhK取此为改进公式。近似公式为 122211212110,1,2( ,)(2,2)nnnnnnccabyyhKKf x yKf xhyhK取此为常用的二阶公式，称为中点公式 三

23、、三阶龙格库塔方法三、三阶龙格库塔方法1123121312(4)6(,)(,)22(,2)nnnnnnnnRKhyyKKKKf xyhhKf xyKKf xh yhKhK常用的三阶公式为：四、四阶龙格库塔方法四、四阶龙格库塔方法112341213243 (22)6(,)(,)22(,)22(,)nnnnnnnnnnRKhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK常用的四阶公式为： 0.2,-83 ;(0)2.0.44 hR Kyyyy例：设取步长写出用经典（标准的）四阶方法求解初值问题 的计算公式，计算的近似值，小数点后至少保留 位。12,8 - 3 , =0.2

24、, 0.2,0.4fx yyhyyyy解：112341122343(22);6,83;,5.62.1;22,6.322.37;22,4.2081.578.nnnnnnnnnnnnnnhyyKKKKKf xyyKhKfxyyKhKfxyyKf xh yKy由经典的四阶龙格库塔公式得 10121.20160.549402,0.22.30040.42.4654nnyyyyyyyy由于11,2,3,4,454652pRKppRKpppRKRKTaylor两点说明：）当时， 公式的最高阶数恰好是 当时， 公式的最高阶数不是，如时仍为 ，时 公式的最高阶数为 。） 方法的导出基于展开，故要求所求问题的解具

25、有较高的光滑度。 RKEulerRKRKEuler当解充分光滑时，四阶 方法确实优于改进法。对一般实际问题，四阶方法一般可达到精度要求。如果解的光滑性差，则用四阶 方法解的效果不如改进法。五、变步长的龙格五、变步长的龙格库塔库塔方法方法( )1( )51115(2)1,(),2,2nhnhnnnnhnxhyy xychhxxhyc 以经典四阶龙格库塔公式为例。从节点 出发，以 为步长求一近似值将步长折半，即取 为步长从 跨两步到，求一近似值每跨一步的截断误差是5(2)11(2)11()11(2)(2)()1111()2,2()1.16()1().15hnnhnnhnnhhhnnnnhy xyc

26、y xyy xyy xyyy因此有由上两式 4 4、微分方程数值解的稳定、微分方程数值解的稳定性性*111*11 nnnnnyyyy稳定性分析，对计算误差：其中是的近似计算值，误差积累会淹没真值？ = , 0 12 nn kn knyyhykh定义：一种数值方法求解 试验方程其中是复常数。对给定的步长，若计算误差在计算， ， 时不产生增大的误差，即，称对与 这种方法是绝对稳定的。 hhh对 ， 的允许范围内是绝对稳定的，则称的全体为绝对稳定区域。Euler法的绝对稳定区域法的绝对稳定区域1*1 nnnnnyyEuleryyhyyhy 的算法：计算值 1 nnnh误差方程：11 11 nnhh

27、从而当是绝对稳定区域 Im(h)-2 -1 0 Re(h)hA绝对稳定区域越大， 可选大些，方法适应性越强。如果整个左半平面是绝对稳定区域称稳定的。向后向后Euler Euler 法的稳定性法的稳定性11 nnnyyEuleryyhy 对用向后法 ： 11 nnnh误差公式：仍受限制。要求稳定的。但收敛法是因此向后 , 10 hhLAEuler1 ()0 1nenRh只 要则11/ 211 111nnhhh122122211112ehhRhh梯形公式的稳定性梯形公式的稳定性11 2nnnyyhnyyyy对用梯形公式）（1()0 1 A- nenRh当时，梯形公式是稳定的。11112212122

28、2 2211()411()4nnnnnnhhheehhRhhRR-K方法的绝对稳定区域方法的绝对稳定区域2121233223443( , ) 11 ,22111 ()22411 ()24nnnnnnnyf x yyRKKh y KhyKyhhKhyKyhhhKhyKyhhhh将代入公式：112342341 22611112624nnnyyKKKKyhhhh234111112624nnhhhh则234111 112624hhhh绝对稳定区域： 2 1-3 -2 -1 0 -1 -211-r1R-K,nnnnnnyyyyyy单步法在计算时，只用到前一步的信息 。为提高精度，需重新计算多个点处的函数

29、值，如方法，计算量较大。如何通过较多地利用前面的已知信息，如 ，来构造高精度的算法计算。基本思想基本思想4.4.线性多步法线性多步法11110111 (,), (,) ,( ,)(1, ,)00Taylornnn rnnn rn rrrnin iin iiiiikkkyyyf xyf xyyyhfff x yknnn r 多步法中最常用的是线性多步法，它的计算公式中只出现， ，及的一次项，其一般形式为 其中均为常数，。若，显式；，隐式。构造线性多步公式常用展开和数值积分方法。线性多步公式的导出1(),nnniiTaylorxTaylory xxTaylor 利用展开导出的基本方法是：将线性多步

30、公式在 处进行展开，然后与在处的展开式相比较，要求它们前面的项重合，由此确定参数。( )( )2( )1() (1,2,),( )( )()()2 ()()!kknnnnnnnnpppnnnyyxky xxTayloryy xyy xxxxyxxO xxp记则在 处的展开为1011110111( ) ()nnnnnnry xyyyhfff以为例：设初值问题的解充分光滑，待定的两步公式为231(4)(5)45(6)21111( ),( )( ,) (),()2!3! ()4!5!(,)()2! iiiiinnnnnnnnnnnnnnnnyy xy xf x yinyyyy xhyy hhhyyh

31、hO hyff xyy xyy hh假设前 步计算结果都是准确的，即则有(4)(5)34(5) ()3!4!nnyyhhO h1111(4)(5)234(5)(,)(,) () ()2!3!4!nnnnnnnnnnnnnff xyyff xyy xyyyyy hhhhO h1(5)2561()()() 2!5!1 nnnnnnnpy xxTayloryyy xyy hhhO hp为使上式有 阶精度，只须使其与在 处的展开式的前项重合。211011101113(4)4111111(5)56111()()()2()()6222466()( )1202424nnnnnnnyyy hy hy hy h

32、y hO h 将以上各公式代入并整理，得010101111111111111122111162261111246624aaaaaa5,5,11,2,3,4)iippp 个参数只须 个条件。由推导知，如果选取参数，使其满足前个方程（，则近似公式为 阶公式。11011111,0,02 ()2nnnnhyyff0如满足方程组前三个方程，故公式此为二阶公式。0111011115(5)6110,1,343 (4)31 ( )90nnnnnnnhyyfffRh yOh又如：解上面方程组得,相应的线性二步四阶公式（Simpson公式)为其截断误差为由此可知，线性二步公式至多是四阶公式。1010123()1(

33、 )( )1 (1,2,)nnnriirrkkiiiirxTaylory xxTaylorikikp一般地，线性多步公式中有个待定参数，如令其右端在 处的展开式与在处的展开式的前p+1项系数对应相等，可得方程组1111(1)21 1( )(1)! (1)( ) ()prpniirpppiniphRippiyO h其解所对应的公式具有 阶精度，局部截断误差为显然，线性多步公式至多可达到2r+2阶精度。二、常用的线性多步公式二、常用的线性多步公式1231010100123 (Adams)r=30,1()()1 (1,2,3 4)5559379=1,24242424riirrkkiiiiikik （

34、一）阿达姆斯公式取，并令由方程组，可解得，1123153354(5)61115(5)6(5559379)24=0Adams1( )5( )()5!251()720nnnnnnniiniinhyyffffhRiiyO hh yO h相应的线性多步公式为因，此式称为显式公式，是四阶公式.局部截断误差为12330101211125(5)610,91951 =1,24242424(9195)24Adams19()720nnnnnnnnhyyffffRh yO h如果令由方程组可解得，相应的线性多步公式为称其为四阶隐式公式，其局部截断误差为利用数值积分方法求线性利用数值积分方法求线性多步公式多步公式11

35、1111()()( , ( )( ),( )( )( ),1nnnnxxnnxxnnn rnnn rry xy xf x y x dxF x dxxxxxxxF xrxF xr 基本思想是首先将初值问题化成等价的积分形式用过节点或的的 次插值多项式代替求积分即得阶的线性多步公式。123330123303,( ) ( )( ) ()()()()()( )()() (0,1,2,3)nnnnin iinnnnin in injjj irx xxxF xL xl x F xxxxxxxxxl xxxxxi例如时，过节点的三次插值多项式为其中1111131301233231313233()()( )(

36、 ) ()()()()()6()()()()2()()()()2()()nnnnnnnnnnxxnnin ixxixnnnnxxnnnnxxnnnnxnny xy xL x dxl x dx F xxxxxxxF xdxhxxxxxxF xdxhxxxxxxF xdxhxxxF x1123123)()655 ()59 ()37 ()9 ()24nnxnnxnnnnxxxdxhhF xF xF xF x1111233,(), (),(,)()(, () (,1,2,3),(5559379)24,nnnnkkkkkknnnnnnnnyyy xy xfxyF xf xy xkn nnnhyyffff

37、AdamsxxAdams对上式用代替用代替则得这就是四阶显式公式。由于积分区间在插值区间外面，又称为四阶外插公式。111(4)310(5)3031(5)35(5)10()()4!() ()4!,),( )251()( )4!720nnnnnnxxnnjxjxxnjxjnnxnnjxjFRxxdxyxxdxxxyRxxdxh y由插值余项公式可得其局部截断误差为由积分中值定理，存在（使得11223111231,( ) ( )( ) ()()()()()( )()() (1,0,1,2)nnnnin iinnnnin in in jjj ixx xxF xL xl x F xxxxxxxxxl x

38、xxxxi 同样，如果过节点的三次插值多项式为其中代( )F x替求积分，11125(5)12121 (9195)2419( )720,nnnnnnnnnnnnAdamshyyffffRh yxxxxAdamsAdams 即得四阶隐式公式其局部截断误差为由于积分区间在插值区间内，故隐式公式又称为内插公式一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程 我们已介绍了一阶常微分方程初值问题的各种我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方种解法，而高阶常微分方程可转化

39、为一阶常微分方程组来求解。程组来求解。 一阶常微分方程组一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组的初值问题 0000( , , ), ()( , , ), ()yf x y z y xyzg x y z z xz（5.15.1） 可以把单个方程可以把单个方程 中的中的f 和和y看作向量看作向量来处理，这样就可把前面介绍的各种差分算法推广来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。到求一阶方程组初值问题中来。 ( , )yf x y 设设 为节点上的近似解，为节点上的近似解，则有改进的则有改进的EulerEuler格式为格式为 0(1,2,3

40、,);,iiixxih iy z1( ,)iiiiiyyhf x y z1( ,)iiiiizzhg x y z预报：预报：1111( ,)(,)2iiiiiiiihyyf x y zf xyz1111( ,)(,)2iiiiiiiihzzg x y zg xyz校正：校正： （7.327.32） 又，相应的四阶龙格又，相应的四阶龙格库塔格式经典格式为库塔格式经典格式为 1123411234(22)6(22)6iiiihyyKKKKhzzLLLL（7.337.33） 112111221112312223122241334133(,)(,)(,)22(,)22(,)22(,)22(,)(,)ii

41、iiiiiiiiiiiiiiiiiiiiiiKfxyzLg xyzhhKfxyKzLhhLg xyKzLhhKfxyKzLhhLg xyKzLKfxyhKzhLLg xyhKzhL式中式中 （7.347.34） 把节点把节点xi上的上的yi和和zi值代入式值代入式(7.34), 依次算出依次算出 , 然后把它们代入式然后把它们代入式(7.33), 算出节点算出节点xi+1上的上的yi+1 和和zi+1值。值。 对于具有三个或三个以上方程的方程组的初值问题对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理也可用类似方法处理,只是更复杂一些而已。此外只是更复杂一些而已。此外,多步多步

42、方法也同样可以应用于求解方程组初值问题。方法也同样可以应用于求解方程组初值问题。 11223344,K L KL KL KL例例 用改进的用改进的EulerEuler法求解初值问题法求解初值问题 (0)1(0)2yxyzyxyzzz 00.2x取步长取步长h=0.1h=0.1，保留六位小数。，保留六位小数。 解解: : 改进的改进的EulerEuler法公式为法公式为),(1iiiiizyxhfyy),(1iiiiizyxhgzz预报：预报： ),(),(21111iiiiiiiizyxfzyxfhyy),(),(21111iiiiiiiizyxgzyxghzz校正：校正： 将将 及及h=0.

43、1h=0.1代入上式代入上式, ,得得 zyxzyxgzxyzyxfiiiiii),(,),(110.1()0.1iiiiiiiiiiyyx yzxyzzz11111110.05 ()()0.05iiiiiiiiiiiiiiiiyyx yzxyzxyxyzzzz由初值由初值 , ,计算得计算得 00(0)1,(0)2yyzz110.8000002.050000yz11(0.1)0.801500(0.1)2.046951yyzz220.6048202.090992yz22(0.2)0.604659(0.2)2.088216yyzz高阶方程组高阶方程组 高阶微分方程高阶微分方程(或方程组或方程组)

44、的初值问题的初值问题,原则上都原则上都可以归结为一阶方程组来求解。例如可以归结为一阶方程组来求解。例如,有二阶微分方有二阶微分方程的初值问题程的初值问题 0000( , ,)(),()yf x y yy xyy xy在引入新的变量在引入新的变量 后后, ,即化为一阶方程组初值问题即化为一阶方程组初值问题: :zy0000( , , ), (), ()zf x y zyz y xy z xy 式式7.367.36为一个一阶方程组的初值问题，对此可为一个一阶方程组的初值问题，对此可用用7.7.17.7.1中介绍的方法来求解。例如应用四阶龙格中介绍的方法来求解。例如应用四阶龙格- -库塔公式得库塔公

45、式得 1123411234(22)6(22)6iiiihyyKKKKhzzLLLL1121211123231222434133(,)2(,)222(,)22(,)iiiiiiiiiiiiiiiiKzLf xy zhKzLhhLf xyKzLhKzLhhLf xyKzLKzhLLf xyhKzhL消去消去 ，上式简化为：，上式简化为： (1,2,3,4)iK i 2112311234()6(22)6iiiiihyyhzLLLhzzLLLL1211223112224123(,)(,)22(,)242(,)2iiiiiiiiiiiiiiiLf xy zhhLf xyzzLhhhLf xyzL zLhLf xyhzLzhL上述方法同样可以用来处理三阶或更高阶的微分方上述方法同样可以用来处理三阶或更高阶的微分方程或方程组