chap常微分方程数值解PPT课件

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

2、令a= x0 x1 xn =b,其中hk=xk+1-xk , 如是等距节点h=(b-a)/n , h称为步长。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk上y(xk)的近似值，用yk表示，即yky(xk)，这样y0 , y1 ,.,yn称为微分方程的数值解。 第4页/共120页主要问题v如何将微分方程离散化，并建立求其数值解的递推公式；v递推公式的局部截断误差，数值解与精确解的误差估计；v递推公式的稳定性与收敛性。第5页/共120页v用差商代替微商用差商代替微商v数值积分数值积分vTaylor展开展开微分方程离散化常用方法微分方程离散化常用方法第6页/

3、共120页第7页/共120页第8页/共120页第9页/共120页第10页/共120页第11页/共120页1 解常微分方程初值问题的Euler方法Euler方法Euler方法的误差分析第12页/共120页v向前向前EulerEuler公式公式（EulerEuler折线法或显格式）折线法或显格式）v向后Euler公式（后退Euler公式）v梯形公式（改进的Euler公式）vEuler预估校正格式一、Euler方法第13页/共120页1、向前Euler公式第14页/共120页第15页/共120页用分段的折线逼近逼近函数第16页/共120页2、向后（后退的）Euler 方法第17页/共120页第18页

4、/共120页3、梯形公式第19页/共120页第20页/共120页4、改进的尤拉公式梯形公式虽然提高了精度，但使算法复杂。而在实际计算中只迭代一次，这样建立的预测校正系统称作改进的尤拉公式。第21页/共120页1111 ( ,); ( ,)(,),2nnnnnnnnnnyyhf x yhyyf x yf xy预测校正11(,);(,);() / 2.pnnncnnpnpcyyhfxyyyhfxyyyy第22页/共120页二、二、Euler方法的误差分析方法的误差分析第23页/共120页第24页/共120页第25页/共120页第26页/共120页2） 总体方法误差总体方法误差第27页/共120页第

5、28页/共120页第29页/共120页第30页/共120页总体截断误差与局部截断误差的关系是：第31页/共120页误差分析表误差分析表EulerEuler方法方法局部截局部截断误差断误差总体截总体截断误差断误差迭代收敛迭代收敛条件条件向前向前EulerEuler方法方法O(h2)O(h)向后向后EulerEuler方法方法O(h2)O(h)0hL1梯形公式梯形公式O(h3)O(h2)0hL2(L为为Lip常数）常数）第32页/共120页向后Euler 方法收敛条件与截断误差第33页/共120页梯形公式的收敛性第34页/共120页2(01);(0)1.xyyxyy例：用尤拉公式和改进的尤拉公式解

6、初值问题10.12().nnnnnhxEuleryyh yy解：取步长，公式为： 第35页/共120页112();2();1().2npnnnncnppnpcxyyh yyxyyh yyyyy改进的尤拉公式为：第36页/共120页20.2(00.6); (0)1.hyyxyxy 例1：取步长，用欧拉法解初值问题122 , 0.2 0.80.2nnnnnnnnnnnEuleryyhfxyyyx yyx y解：格式为：第37页/共120页0122 1 0.2 0.8, 0.4 0.6144 0.6 0.461321yyyyyyy由计 算 得第38页/共120页0.283(12); (1)2.5hy

7、yxy例2：取步长，用梯形解初值问题小数点后至少保留 位。11111 ,20.2 83832nnnnnnnnnnhyyfxyfxyyyyy解：梯形公式为： 第39页/共120页 1012345716 131312, 1.22.307691.42.473371.62.562581.82.610622.02.63649nnyyyyyyyyyyyyyy故 由计 算 得 第40页/共120页0.23(12); (1)2.5hyxyxy例2 ：取步长，用梯形解初值问题小数点后至少保留 位。111111 ,2 0.1 33nnnnnnnnnnnnhyyfxyfxyyyx yxy解：梯形公式为：第41页/共

8、120页 0111110011234511111110.30.1 332,2.6*,&,%,?nnnnkknnnnnnkyyx yyyx yxyyyyyyyyyy迭代格式：第42页/共120页 021234522222,*,&,%,yyyyyy第43页/共120页20.2sin0 (1)1.1.21.45hyyyxyyy例3：取步长，用欧拉预校方法解初值问题计算及的近似值，小数点后至少保留 位。第44页/共120页111121212111,2 0.2sin0.1sin + sinnnnnnnnnnnnnnnnnnnnnnnnyyhfxyhyyfxyfxyyyyyxyyyyxyy

9、x解：欧拉预校格式为：第45页/共120页 0112211,0.63171 1.20.7154880.476961.40.52611yyyyyyyy故 由计 算 得 第46页/共120页50 (0)1.0.10.210yyyhy例4：写出用反复迭代的欧拉预校法解初值问题的计算公式，并取步长，计算，要求迭代误差不超过。第47页/共120页 01111101111,2 0,1,2, 0,1,2, 0.90.950.05nnnnkknnnnnnnnkknnnyyhfxyhyyfxyfxyknfx yyyyyyy 解：欧拉预校的迭代格式为：取第48页/共120页 0012111341143751141

10、101, 0.9,0.905,0.90475, 0.9047625,0.9047618756.251010, 0.10.904761875yyyyyyyyyyyy故 由计 算 得由于 是 取第49页/共120页 0012232452254652252201, 0.814286,0.818809, 0.818583, 0.818595,0.8185941010, 0.20.818594yyyyyyyyyyyy故 由计 算 得由于 是 取第50页/共120页2.2.龙格龙格库塔方法库塔方法基本思想基本思想二阶二阶R-K方法方法三阶三阶R-K方法方法四阶四阶R-K方法方法变步长变步长R-K方法方法第

11、51页/共120页1112() ()()()()()()2!( )( , ), ( )( , )( , ) ( , ),nnnppnnnnxypTaylory xyy xhhy xhy xyxyxPy xf x yyxfx yfx y f x y若用 阶多项式近似函数有：其中。但由于公式中各阶偏导数计算复杂，不实用。一、基本思想第52页/共120页(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1); ;2,3,jjjjyffffyffxyffyffxyffyffjxy一般地有第53页/共120页111112121(,)11()22 (,)(,)nnnnnnnnnnEulerEu

12、lerEuleryyhKKf xyEuleryyhKKKf xyKf xh yhK如果将公式与改进公式写成下列形式：公式改进公式第54页/共120页11 ( , )()( , )( , )nnf x yy xyf x yf x y以上两组公式都使用函数在某些点上的值的线性组合来计算的近似值。Euler公式：每步计算一次的值，为一阶方法。改进Euler公式：需计算两次的值，二阶方法。第55页/共120页( , )(,)( )-nnnf x yxyTaylory xxTaylorR K于是可考虑用函数在若干点上的函数值的线性组合来构造近似公式，构造是要求近似公式在处的展开式与解在 处的展开式的前面

13、几项重合，从而使近似公式达到所需要的阶数。即避免求偏导，又提高了方法的精度，此为方法的基本思想。第56页/共120页11111-(,)(,) (2,3,),pnniiinniininijjjiijiRKyyhc KKf xyKf xa h yhb Kipabc一般地，方法设近似公式为其中 ，都是参数.第57页/共120页,(,)( )iijinnnab cxyTaylory xxTaylor确定参数 ，的原则是:使近似公式在处的展开式与在 处的展开式的前面项尽可能多地重合。第58页/共120页二、二阶龙格库塔方法111221222112()(,)(,)nnnnnnpyyh c Kc KKf x

14、yKf xa h yhb K当时，近似公式为 第59页/共120页112221123221( ,)( ,)(, ( ,)( ,) ( ,) ( ,)( ,) ( ,)( )nnnnnnnnnnnnnnnxnnynnnnnx 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 h第60页/共120页123123()()

15、()()()()2(,)(,)2 (,) (,)()nnnnnnnnnxnnynnnny xxTaylorhy xy xhy xy xO hhyf xy hfxyfxyf xyO h在 处的展开式为第61页/共120页122 232 211 1/2 1/2 ( ),cccacbOh有 无 穷 多 组 解 ， 每 一 组 解 得 一近 似 公 式 ， 局 部 截 断 误 差 均 为这 些 方 法 统 称 二 阶 方 法 。43()()O hO h可以证明，无论这四个参数如何选择，都不能使局部截断误差达到，也即在计算两次函数值的情况下，局部截断误差的阶最高为。第62页/共120页122211121

16、211,1,2()/2(,)(,)nnnnnnccabEuleryyh KKKf xyKf xh yhK取此为改进公式。近似公式为 第63页/共120页122211212110,1,2( ,)(2,2)nnnnnnccabyyhKKf x yKf xhyhK取此为常用的二阶公式，称为中点公式 第64页/共120页三、三阶龙格库塔方法1123121312(4)6(,)(,)22(,2)nnnnnnnnRKhyyKKKKf xyhhKf xyKKf xh yhKhK常用的三阶公式为：第65页/共120页四、四阶龙格库塔方法112341213243 (22)6(,)(,)22(,)22(,)nnnn

17、nnnnnnRKhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK常用的四阶公式为： 第66页/共120页0.2,-83 ;(0)2.0.44 hR Kyyyy例：设取步长写出用经典（标准的）四阶方法求解初值问题 的计算公式，计算的近似值，小数点后至少保留 位。12,8 - 3 , =0.2, 0.2,0.4fx yyhyyyy解：第67页/共120页112341122343(22);6,83;,5.62.1;22,6.322.37;22,4.2081.578.nnnnnnnnnnnnnnhyyKKKKKf xyyKhKfxyyKhKfxyyKf xh yKy由经典的四

18、阶龙格库塔公式得第68页/共120页 10121.20160.549402,0.22.30040.42.4654nnyyyyyyyy由于第69页/共120页11,2,3,4,454652pR KppR KpppR KR KTaylor两点说明：）当时， 公式的最高阶数恰好是 当时， 公式的最高阶数不是，如时仍为 ，时 公式的最高阶数为。） 方法的导出基于展开，故要求所求问题的解具有较高的光滑度。第70页/共120页 RKEulerRKRKEuler当解充分光滑时，四阶 方法确实优于改进法。对一般实际问题，四阶方法一般可达到精度要求。如果解的光滑性差，则用四阶 方法解的效果不如改进法。第71页/

19、共120页五、变步长的龙格库塔方法( )1( )51115(2)1,(),2,2nhnhnnnnhnxhyy xychhxxhyc 以经典四阶龙格库塔公式为例。从节点 出发，以 为步长求一近似值将步长折半，即取 为步长从 跨两步到，求一近似值每跨一步的截断误差是第72页/共120页5(2)11(2)11()11(2)(2)()1111()2,2()1.16()1().15hnnhnnhnnhhhnnnnhy xycy xyy xyy xyyy因此有由上两式 第73页/共120页 4、微分方程数值解的稳定性第74页/共120页第75页/共120页Euler法的绝对稳定区域第76页/共120页 I

20、m(h)-2 -1 0 Re(h)第77页/共120页向后向后Euler Euler 法的稳定性法的稳定性第78页/共120页第79页/共120页梯形公式的稳定性第80页/共120页第81页/共120页R-K方法的绝对稳定区域第82页/共120页第83页/共120页 2 1-3 -2 -1 0 -1 -2第84页/共120页11-r1R-K,nnnnnnyyyyyy单步法在计算时，只用到前一步的信息 。为提高精度，需重新计算多个点处的函数值，如方法，计算量较大。如何通过较多地利用前面的已知信息，如 ，来构造高精度的算法计算。基本思想4.线性多步法第85页/共120页11110111 (,),

21、(,) ,( ,)(1, ,)00Taylornnn rnnn rn rrrnin iin iiiiikkkyyyf xyf xyyyhfff x yknnn r 多步法中最常用的是线性多步法，它的计算公式中只出现， ，及的一次项，其一般形式为 其中均为常数，。若，显式；，隐式。构造线性多步公式常用展开和数值积分方法。第86页/共120页线性多步公式的导出1(),nnniiTaylorxTaylory xxTaylor 利用展开导出的基本方法是：将线性多步公式在 处进行展开，然后与在处的展开式相比较，要求它们前面的项重合，由此确定参数。第87页/共120页( )( )2( )1() (1,2,

22、),( )( )()()2 ()()!kknnnnnnnnpppnnnyyxky xxTayloryy xyy xxxxyxxO xxp记则在 处的展开为1011110111( ) ()nnnnnnry xyyyhfff以为例：设初值问题的解充分光滑，待定的两步公式为第88页/共120页231(4)(5)45(6)21111(),()( ,) (),()2!3! ()4!5!(,)()2! iiiiinnnnnnnnnnnnnnnnyy xy xf x yinyyyy xhyy hhhyyhhO hyff xyy xyy hh假设前 步计算结果都是准确的，即则有(4)(5)34(5) ()3!

23、4!nnyyhhO h第89页/共120页1111(4)(5)234(5)(,)(,) () ()2!3!4!nnnnnnnnnnnnnff xyyff xyy xyyyyy hhhhO h第90页/共120页1(5)2561()()() 2!5!1 nnnnnnnpy xxTayloryyy xyy hhhO hp为使上式有 阶精度，只须使其与在 处的展开式的前项重合。211011101113(4)4111111(5)56111()()()2()()6222466()( )1202424nnnnnnnyyy hy hy hy hy hO h 将以上各公式代入并整理，得第91页/共120页01

24、0101111111111111122111162261111246624aaaaaa5,5,11,2,3,4)iippp 个参数只须 个条件。由推导知，如果选取参数，使其满足前个方程（，则近似公式为 阶公式。第92页/共120页11011111,0,02 ()2nnnnhyyff0如满足方程组前三个方程，故公式此为二阶公式。第93页/共120页0111011115(5)6110,1,343 (4)31 ( )90nnnnnnnhyyfffRh yOh又如：解上面方程组得,相应的线性二步四阶公式（Simpson公式)为其截断误差为由此可知，线性二步公式至多是四阶公式。第94页/共120页101

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

26、1,24242424riirrkkiiiiikik （一）阿达姆斯公式取，并令由方程组，可解得，第98页/共120页1123153354(5)61115(5)6(5559379)24=0Adams1( )5( )()5!251()720nnnnnnniiniinhyyffffhRiiyO hh yO h相应的线性多步公式为因，此式称为显式公式，是四阶公式.局部截断误差为第99页/共120页12330101211125(5)610,91951 =1,24242424(9195)24Adams19()720nnnnnnnnhyyffffRh yO h如果令由方程组可解得，相应的线性多步公式为称其为

27、四阶隐式公式，其局部截断误差为第100页/共120页利用数值积分方法求线性多步公式111111()()( , ( )( ),( )( )( ),1nnnnxxnnxxnnn rnnn rry xy xf x y x dxF x dxxxxxxxF xrxF xr 基本思想是首先将初值问题化成等价的积分形式用过节点或的的 次插值多项式代替求积分即得阶的线性多步公式。第101页/共120页123330123303,( ) ( )( ) ()()()()()( )()() (0,1,2,3)nnnnin iinnnnin in injjj irx xxxF xL xl x F xxxxxxxxxl

28、xxxxxi例如时，过节点的三次插值多项式为其中第102页/共120页1111131301233231313233()()( )( ) ()()()()()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 x第103页/共120页1111233,(), (),(,)()(,

29、 () (,1,2,3),(5559379)24,nnnnkkkkkknnnnnnnnyyy xy xfxyF xf xy xkn nnnhyyffffAdamsxxAdams对上式用代替用代替则得这就是四阶显式公式。由于积分区间在插值区间外面，又称为四阶外插公式。第104页/共120页111(4)310(5)3031(5)35(5)10()()4!() ()4!,),( )251()( )4!720nnnnnnxxnnjxjxxnjxjnnxnnjxjFRxxdxyxxdxxxyRxxdxh y由插值余项公式可得其局部截断误差为由积分中值定理，存在（使得第105页/共120页11223111

30、231,( ) ( )( ) ()()()()()( )()() (1,0,1,2)nnnnin iinnnnin in in jjj ixx xxF xL xl x F xxxxxxxxxl xxxxxi 同样，如果过节点的三次插值多项式为其中代( )F x替求积分，第106页/共120页11125(5)12121 (9195)2419( )720,nnnnnnnnnnnnAdamshyyffffRh yxxxxAdamsAdams 即得四阶隐式公式其局部截断误差为由于积分区间在插值区间内，故隐式公式又称为内插公式第107页/共120页一阶常微分方程组与高阶方程 我们已介绍了一阶常微分方程初

31、值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。一阶常微分方程组对于一阶常微分方程组的初值问题 0000( , , ), ()( , , ), ()yf x y z y xyzg x y z z xz（） 可以把单个方程 中的f 和y看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。 ( , )yf x y 第108页/共120页设 为节点上的近似解，则有改进的Euler格式为 0(1,2,3,);,iiixxih iy z1( ,)iiiiiyyhf x y z1( ,)iiiiizzhg x y

32、z预报：1111( ,)(,)2iiiiiiiihyyf x y zf xyz1111( ,)(,)2iiiiiiiihzzg x y zg xyz校正： （） 又，相应的四阶龙格库塔格式（经典格式）为 1123411234(22)6(22)6iiiihyyKKKKhzzLLLL（） 第109页/共120页112111221112312223122241334133(,)(,)(,)22(,)22(,)22(,)22(,)(,)iiiiiiiiiiiiiiiiiiiiiiiiKfxyzLg xyzhhKfxyKzLhhLg xyKzLhhKfxyKzLhhLg xyKzLKfxyhKzhLLg

33、 xyhKzhL式中 （） 第110页/共120页 把节点xi上的yi和zi值代入式(7.34), 依次算出 , 然后把它们代入式(7.33), 算出节点xi+1上的yi+1 和zi+1值。 对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。 11223344,K L KL KL KL例 用改进的Euler法求解初值问题 (0)1(0)2yxyzyxyzzz 00.2x取步长，保留六位小数。 解: 改进的Euler法公式为),(1iiiiizyxhfyy),(1iiiiizyxhgzz预报： ),(),(21

34、111iiiiiiiizyxfzyxfhyy),(),(21111iiiiiiiizyxgzyxghzz校正： 将 及代入上式,得 zyxzyxgzxyzyxfiiiiii),(,),(第111页/共120页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第112页/共120页高阶方程组 高阶微分方程(或方程组)的初值问题,原则上都可以归结为一阶方程组来求解。例如,有二阶微分方程的初值问题 0000( , ,)(),()yf x y yy xyy xy在引入新的变量 后,即化为一阶方程组初值问题:zy0000( , , ), (), ()zf x y zyz y xy z xy 式（）为一个一阶方程组的初值问题，对此可用中介绍的方法来求解。例如应用四阶龙格-库塔