9常微分方程数值解59051

1、/* numerical methods for ordinary differential equations */ 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 /* initial-value problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b r1 上连续，且关于上连续，且关于 y 满足满足 lipschitz 条条件件，即存在与，即存在与 x, y 无关的常数无关的常数 l 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，则上述ivp存存在唯一解在唯一解。| ),(

2、),(|2121yylyxfyxf 要计算出解函数要计算出解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距 为步长，通常采用为步长，通常采用等距节点等距节点，即取即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii第第9章章 常微分方程数值解常微分方程数值解1 引言引言2 欧拉方法欧拉方法/* eulers method */ 欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy

3、1y记为记为)1,., 0(),(1 niyxfhyyiiii定义定义在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 ri = y(xi+1) yi+1 称为称为局部截断误差局部截断误差 /* local truncation error */。定义定义若某算法的局部截断误差为若某算法的局部截断误差为o(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()()(32112iiiihiiiiiyxhfyhoxyxyhxyyxyr )()(322ho

4、xyih 欧拉法具有欧拉法具有 1 阶精度。阶精度。ri 的的主项主项/* leading term */亦称为亦称为欧拉折线法欧拉折线法 /* eulers polygonal arc method*/ 例例1：解：解：h=0.2 , xi=1+ih， 精确解为：精确解为：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 dy2 y2dxxy 11 可

5、以看出误差随着计算在积累。可以看出误差随着计算在积累。xky(xk)ykek1.2- 0.96-10.041.4-0.84-0.93330.09331.6-0.64-0.80.161.8-0.36-0.60.242.00-0.33330.33332.20.4400.44步长步长h=0.21 eulers method 隐式欧拉法隐式欧拉法 /* implicit euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时

6、出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(iiiyxyr)()(322hoxyih 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 欧拉公式的改进：欧拉公式的改进：1 eulers method/* trapezoid formula */ 显显/隐式两种算法的隐式两种算

7、法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法。迭代法。)()(311hoyxyriii 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyii

8、ii假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hoyxyriii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。 梯形公式梯形公式方方 法法 1 eulers method显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度

9、低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 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./* modified eulers method */step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiii

10、yxfhyy step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。) 1,., 0(),(,),(211niyxfhyxf

11、yxfhyyiiiiiiii1 eulers method 改进欧拉法改进欧拉法例例2：解：梯形公式为：解：梯形公式为：, ( )dy2xyy 01dxykk 1k 1kkk 1kk 1h2x2xyyyy2yy ()()2kk 1k 1kk 1khhhxy1y1yhx022y 可以预测校正方法来求：可以预测校正方法来求：()()()kpkkkk 1ckppk 1pc2xyyh yy2xyyh yy1yyy2 解解yk+1出来比较困难，遇到的是一个二次方程，出来比较困难，遇到的是一个二次方程，即：即：),(,),(211kkkkkkkkyxfhyxfyxfhyyxkyky(xk)ek0.11.0

12、9596911.09544510.0004640.21.18409661.18321600.0008810.31.26620141.26491110.0012900.41.34336021.34164080.0017190.51.41640191.41421360.0021880.61.47595561.48323970.0027160.71.55251411.54919330.0033210.81.61647481.61245150.0040230.91.67816641.67332010.0047461.01.73786741.73205080.005817),(,),(211kkkkkk

13、kkyxfhyxfyxfhyy/* runge-kutta method */建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hkyhxfkyxfkkkhyyiiiiii 斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？步长一定是

14、一个步长一定是一个h 吗？吗？3 龙格龙格 - 库塔法库塔法2 runge-kutta method首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hoyxyriii step 1: 将将 k2 在在 ( xi , yi ) 点作点作 taylor 展开展开)(),(),(),(),(2112hoyxfphkyxphfyxfphkyphxfkiiyiixiiii )()()(2hoxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(1212

15、2111phkyphxfkyxfkkkhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx step 2: 将将 k2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hoxyphxyhyhoxyphxyxyhyyiiiiiiii 2 runge-kutta methodstep 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hoxyphxyhyyiiii )()(2)()()(321hoxyhxyhxyxyiii

16、i 要求要求 ，则必须有：，则必须有：)()(311hoyxyriii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些

17、系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 2 runge-kutta method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihkhkhkyihxfkhkhkyhxfkhkyhxfkyxfkkkkhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* classical runge-kutta method */ ：),(),(),(),()22(34222312221432161hkyhxfkkyxfkkyxfkyxfkkkkkyyiihihihihiiihii 例例3：0021

18、737. 222222,22000202003khyhxhxkhykhyhxfk007603. 2,3000303004hkyhxhxhkyhkyhxfk 2002714. 1 326432101kkkkhyy 11yyxxydxdy 1ln2xxxy其精确解为：其精确解为：用经典龙格库塔方法用经典龙格库塔方法，h=0.11 . 0, 1, 100hxy(1) 求求y1，此时，此时2,0000001yxxyyxfk0021645. 222222,21000101002khyhxhxkhykhyhxfk解：解： 2002711. 11ln2111xxxy以下计算用表格列出：以下计算用表格列出：

19、kxkykxyke7103710571061.11.20027141.20027111.21.40181581.40181531.31.60523931.60523871.41.81079361.71079301.52.01856282.01856211.62.22854702.22854631.72.44070382.44070301.82.65496932.65496851.92.87126972.87126892.03.08952793.08952717106710771077108710871087108附：附：三阶龙格库塔方法三阶龙格库塔方法:122211333113221112233(,)(,)(,()()kkkkkkkkkf xykf xa h yhkkf xa h yhkkyyhkkk 共有八个参数：共有八个参数： 1, 2, 3, 1, 2, 21, 31, 32，1、三阶龙格库塔公式格式如下：、三阶龙格库塔公式格式如下：v 将将 f 在在 xk 点点泰勒展开泰勒展开： 1,kkkkf x yy x 2222223222122212111()22xyxxxyyykfa hfhffa h fah ffh f fo h 22233311