常微分方程数值解59382PPT学习教案

1、会计学1 常微分方程数值解常微分方程数值解59382 /* Eulers Method */ 欧拉公式：欧拉公式： x0 x1 向前差商近似导数向前差商近似导数 h xyxy xy )()( )( 01 0 ),()()()( 000001 yxfhyxyhxyxy 1 y 记为记为 )1,., 0(),( 1 niyxfhyy iiii 定义定义在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考步计算是精确的前提下，考 虑的截断误差虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断误差局部截断误差 /* local truncation error

2、*/。 定义定义若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 欧拉法的局部截断误差：欧拉法的局部截断误差： ),()()()()()( 3 211 2 iiii h iiiii yxhfyhOxyxyhxyyxyR )()( 3 2 2 hOxy i h 欧拉法具有欧拉法具有 1 阶精度。阶精度。 Ri 的的主项主项 /* leading term */ 亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 第1页/共19页 解：解：h=0.2 , xi=1+ih， 精确解为：精确解

3、为：y=x2-2x y1= y0+hf(x0, y0)=-1 y2= y1+hf(x1, y1)=- 0.9333 y3= y2+hf(x2, y2)=-0.8 y4= y3+hf(x3, y3)=-0.6 y5= y4+hf(x4, y4)=-0.3333 y6= y5+hf(x5, y5)=0 dy2y 2 dxx y 11 可以看出误差随着计算在积累。可以看出误差随着计算在积累。 xky(xk)ykek 1.2- 0.96-10.04 1.4-0.84 - 0.9333 0.0933 1.6-0.64-0.80.16 1.8-0.36-0.60.24 2.00 - 0.3333 0.33

4、33 2.20.4400.44 步长步长h=0.2 第2页/共19页 1 Eulers Method 隐式欧拉法隐式欧拉法 /* implicit Euler method */ 向后差商近似导数向后差商近似导数 h xyxy xy )()( )( 01 1 x0 x1 )(,()( 1101 xyxfhyxy )1,., 0(),( 111 niyxfhyy iiii 由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故 称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explic

5、it */ 欧拉公式。欧拉公式。 一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差： 11) ( iii yxyR)()( 3 2 2 hOxy i h 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 第3页/共19页 1 Eulers Method /* trapezoid formula */ 显显/隐式两种算法的隐式两种算法的平均平均 )1,., 0(),(),( 2 111 niyxfyxf h yy iiiiii 注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2

6、 阶精度，比欧拉方法有了进阶精度，比欧拉方法有了进 步。但注意到该公式是步。但注意到该公式是隐式隐式公式，计算时不得公式，计算时不得 不用到迭代法。不用到迭代法。 )()( 3 11 hOyxyR iii 中点欧拉公式中点欧拉公式 /* midpoint formula */ 中心差商近似导数中心差商近似导数 h xyxy xy 2 )()( )( 02 1 x0 x2x1 )(,(2)()( 1102 xyxfhxyxy 1,., 1),(2 11 niyxfhyy iiii 假设假设 ，则可以导出，则可以导出 即中点公式具有即中点公式具有 2 阶精度。阶精度。 )(),( 11iiii x

7、yyxyy )()( 3 11 hOyxyR iii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推 过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。 第4页/共19页 方方 法法 1 Eulers Method 显式欧拉显式欧拉 隐式欧拉隐式欧拉 梯形公式梯形公式 中点公式中点公式 简单简单精度低精度低 稳定性最好稳定性最好精度低精度低, 计算量大计算量大 精度提高精度提高计算量大计算量大 精度提高精度提高, 显式显式

8、多一个初值多一个初值, 可能影响精度可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets make it possible. 第5页/共19页 /* modified Eulers method */ Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出 ),( 1iiii yxfhyy Step 2: 再将再将 代入代入隐式隐

9、式梯形公式的右边作梯形公式的右边作校正校正，得到，得到 1 i y ),(),( 2 111 iiiiii yxfyxf h yy 注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。 可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个 单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后。后 面将看到，它的面将看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。 ) 1,., 0(),(,),( 2 11 niyxfhyxfyxf h yy ii

10、iiiiii 1 Eulers Method 第6页/共19页 解：梯形公式为：解：梯形公式为： , ( ) dy2x yy 01 dxy kk 1 k 1kkk 1 kk 1 h2x2x yyyy 2yy ()() 2k k 1 k 1kk 1 k hhhx y1y1yhx0 22y 可以预测校正方法来求：可以预测校正方法来求： () () () k pkk k k 1 ckp p k 1pc 2x yyh y y 2x yyh y y 1 yyy 2 解解yk+1出来比较困难，遇到的是一个二次方程，出来比较困难，遇到的是一个二次方程， 即：即： ),(,),( 2 11kkkkkkkk y

11、xfhyxfyxf h yy 第7页/共19页 xkyky(xk)ek 0.11.09596911.09544510.000464 0.21.18409661.18321600.000881 0.31.26620141.26491110.001290 0.41.34336021.34164080.001719 0.51.41640191.41421360.002188 0.61.47595561.48323970.002716 0.71.55251411.54919330.003321 0.81.61647481.61245150.004023 0.91.67816641.67332010.0

12、04746 1.01.73786741.73205080.005817 ),(,),( 2 11kkkkkkkk yxfhyxfyxf h yy 第8页/共19页 /* Runge-Kutta Method */ 建立高精度的单步递推格式。建立高精度的单步递推格式。 单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜 率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所 能达到的最高精度为能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：

13、),( ),( 2 1 2 1 12 1 211 hKyhxfK yxfK KKhyy ii ii ii 斜率斜率 一定取一定取K1 K2 的的平均值平均值吗？吗？ 步长一定是一个步长一定是一个h 吗吗 ？ 第9页/共19页 2 Runge-Kutta Method 首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶 精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )( ii xyy )()( 3 11 hOyxyR iii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开 )(),(),(),

14、( ),( 2 1 12 hOyxfphKyxphfyxf phKyphxfK iiyiixii ii )()()( 2 hOxyphxy ii 将改进欧拉法推广为：将改进欧拉法推广为： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii ),(),(),( ),(),( ),()( yxfyxfyxf dx dy yxfyxf yxf dx d xy yx yx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()( )()()()( 32 221 2 211 hOxyphxyhy hOxyphxyxyhyy iii iii

15、ii 第10页/共19页 2 Runge-Kutta Method Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较 )()()()( 32 2211 hOxyphxyhyy iiii )()( 2 )()()( 3 2 1 hOxy h xyhxyxy iiii 要求要求 ，则必须有，则必须有 ： )()( 3 11 hOyxyR iii 2 1 ,1 221 p 这里有这里有 个未知个未知 数，数， 个方程。个方程。 3 2 存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库塔格式

16、库塔格式。 2 1 , 1 21 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？ 第11页/共19页 其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为均为 待定系数，确定这些系待定系数，确定这些系 数的步骤与前面相似。数的步骤与前面相似。 2 Runge-Kutta Method ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 mm mmmmi

17、m ii ii ii mmii hKhKhKyihxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ： ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK KyxfK yxfK KKKKyy ii h i h i h i h i ii h ii 第12页/共19页 0021737. 2 2 2 2 2 2 , 2 20 0 0 20 2003 K h y h x h x K h y K

18、 h y h xfK 007603. 2, 30 0 0 30 3004 hKy hx hx hKy hKyhxfK 2002714. 1 32 6 432101 KKKK h yy 11y y x x y dx dy 1ln2xxxy 其精确解为：其精确解为： 用经典龙格库塔方法用经典龙格库塔方法，h=0.1 1 . 0, 1, 1 00 hxy(1) 求求y1，此时，此时 2, 0 0 0 0 001 y x x y yxfK 0021645. 2 2 2 2 2 2 , 2 10 0 0 10 1002 K h y h x h x K h y K h y h xfK 解：解： 20027

19、11. 11ln2 111 xxxy 第13页/共19页 k x k y k xy k e 7 103 7 105 7 106 1.11.200271 4 1.2002711 1.21.401815 8 1.4018153 1.31.605239 3 1.6052387 1.41.810793 6 1.7107930 1.52.018562 8 2.0185621 1.62.228547 0 2.2285463 1.72.440703 8 2.4407030 1.82.654969 3 2.6549685 1.92.871269 7 2.8712689 2.03.089527 9 3.0895

20、271 7 106 7 107 7 107 7 108 7 108 7 108 7 108 第14页/共19页 1 22211 33311322 1112233 (,) (,) (,() () kk kk kk kk Kf xy Kf xa h yhK Kf xa h yhKK yyhKKK 共有八个参数：共有八个参数： 1, 2, 3, 1, 2, 21, 31, 32， 1、三阶龙格库塔公式格式如下：、三阶龙格库塔公式格式如下： v 将将 f 在在 xk 点点泰勒展开泰勒展开： 1 , kkk Kf xyyx 2222223 2221222121 11 () 22 xyxxxyyy Kfa hfhffa h fah ffh f fO h 222 3331132233311322 2 3 311322 1 () 2 1 2 xyxxxy yy KfhfKKhffh fKKh f KKfh 第15页/共19页 将将K1,K2按展开式代入按展开式