CH01常微方程数值解解析

1、.1 .2 主要内容 第第3 3章章 抛物型方程的有限差分方法抛物型方程的有限差分方法 第第1 1章章 常微分方程初值问题数值解法常微分方程初值问题数值解法 第第2 2章章 椭圆型方程的有限差分方法椭圆型方程的有限差分方法 第第4 4章章 双曲型方程的有限差分方法双曲型方程的有限差分方法 第第5 5章章 非线性双曲型守恒律方程的差分方法非线性双曲型守恒律方程的差分方法 第第6 6章章 有限元方法简介有限元方法简介 .3 第第1 1章章 常微分方程初值问题的常微分方程初值问题的 数值解法数值解法 1.11.1 引言引言 1.2 1.2 欧拉法欧拉法( (Euler方法方法) 1.3 1.3 梯形

2、法、隐式格式的迭代计算梯形法、隐式格式的迭代计算 1.4 1.4 一般单步法、一般单步法、Runge-Kutta格式格式 1.5 1.5 线性多步法线性多步法 1.6 1.6 误差的事后估计法、步长的自动选择误差的事后估计法、步长的自动选择 1.7 1.7 高阶常微分方程高阶常微分方程( (组组) )的数值方法的数值方法 .4 1.1 1.1 引言引言 目标在于给出解在一些离散点上的近似值。目标在于给出解在一些离散点上的近似值。 本章研究常微分方程初值问题的主要数值本章研究常微分方程初值问题的主要数值 解法，包括基本方法和基本理论问题。解法，包括基本方法和基本理论问题。 .5 1.2 1.2

3、欧拉法欧拉法( (Euler方法方法) ) 1.2.1 1.2.1 欧拉方法欧拉方法 )(| ),( 0 0 0 xyy Xxxyxf dx dy xx )2 . 1 ( ) 1 . 1 ( 考虑常微分方程初值问题 注：在后面的讨论中，我们总认为这个初值问 题的解存在、唯一且连续依赖于初值条件，即 初值问题(1.1)，(1.2)是适定的。 .6 ),(yxf dx dy )(,()()( 1nnnn xyxhfxyxy 将解的存在区间 N 等分，得到N 个小区间。 任取一个小区间 , 1nn xx，由原方程 1 ),()()( 1 n n x x nn dxyxfxyxy得 ),(yxf在区间

4、 , 1nn xx上，用 在点 nx上的值来 代替),(yxf，得到 N xX xxh nn 0 1 为步长。 其中 .7 图图1.1 0 )(,(xyxf n x 1n xx .8 .)( ! 2 1 )( )()( 2 1 nnnn xyhxhyxyxy )(,()( nnn xyxfxy )(,()()( 1nnnn xyxhfxyxy 在上式中分别用 n y和来代替 1n y)( n xy)( 1n xy和 并由n 的任意性，得到 )( ,.2 , 1 , 0),( 00 1 xyy nyxhfyy nnnn (1.3) 这就是欧拉公式。 欧拉公式亦可由Taylor 展式得到 .9 在

5、上式中分别用 n y和来代替 1n y)( n xy)( 1n xy和 ),( 1nnnn yxhfyy 则得 ),()( 11nnnnnn yxfxxyy 一般而言，并不要求步长相等，则有 (1.4) .10 0 图图1.2 n x y 0 y 3210 xxxx xyy x 几何意义 .11 例例 1.11.1 以 的数值解，并与精确解 为步长，用欧拉法求初值问题1 . 0h 1)0(y yxe dx dy x 比较。 x exxy )2( 2 1 )( 2 .12 (1)、计算格式本身不能准确描述原来的方程 误差的产生：误差的产生： (2)、计算机本身引入的误差 .13 计算机输出的是欧

6、拉方程的近似解 ,而不是 精确解 。因此 n y )( n xy ) ()()( nnnnnn yyxyyxyy (1.5) 可见，为了使计算得到的解 是 的好的近 似，我们要求： n y )( n xy (1)欧拉方法的精确解 是微分方程精确解 的很好近似。 n y)( n xy (2) 是 好的近似。 n y n y 问题(1)称为格式的收敛性问题。 问题(2)称为格式的稳定性问题。 .14 1.2.2 收敛性研究收敛性研究 0hxnhxxn 0 )(xyyn 所谓收敛性问题,就是研究时, 要求 , 。 111 )( nnn yxy整体截断误差整体截断误差 局部截断误差局部截断误差 n n

7、 nn Ryxye1 * 11 )( )(,()(1 * nnn nxyxhfxyy这里 1 ),()(,()()( 11 n n x x nnnnnn yxhfydxxyxfxyyxy(1.6) 即 Th Th .15 定理定理 1.11.1 2 2 1 |MhR n (1.12) | )(|max 0 xyM Xxx 其中h。为步长， XxCxyy,)( 0 2 n R的局部截断误差 ，则欧拉方法 满足 假定 Back .16 其中 R 为局部截断误差的上界。 定理定理 1.21.2 设 f(x,y) 关于 y 满足Lipschitz 条件， L为相应的Lipschitz常数，则欧拉方法的

8、整体截断 满足 n ) 1(| )( 0 )( 00 LxXLxX n e hL R e(1.13) 误差 Th1.4 .17 由定理1.1，1.2，可得 定理定理 1.31.3 设 f(x,y) 关于 y 满足Lipschitz 条件， L为相应的Lipschitz常数， (1.14) XxCxyy,)( 0 2 )( 00 xyy )( n xy 且当h0 ， ，并有估计式 ，则欧拉方法的解一致收敛到初值问 题(1.1)，(1.2)的解 ) 1( 2 | )( 0 )( 00 xXLxXL n e L Mh e n y )( 00 xyy 0 0 ) 1( 2 | )( 0 xXL n e

9、M L h )(|hO n 如果，即 ，由此有 (1.15) 即 .18 欧拉方法的整体截断误差与欧拉方法的整体截断误差与h h 同阶，由同阶，由 的的 表达式可知，表达式可知， ，这说明局部截断误差比整，这说明局部截断误差比整 体截断误差高一阶。体截断误差高一阶。 n R )( 2 hORn 我们称欧拉方法为一阶格式。我们称欧拉方法为一阶格式。 .19 1.2.3 1.2.3 稳定性研究稳定性研究 前已指出欧拉方法的稳定性问题是决定欧拉法在利 用计算机能否得到精确解的关键问题，只有稳定的算法 才可能是有用的算法。 .20 定义1.1 00,z y 如果存在正常数 c 及 ，使对任意初 始值

10、0 h 0 1 ),( z zxhfzz nnnn 与 0 1 ),( y yxhfyy nnnn ，由 nn zy ,计算所得之解满足估计式 0000 ;0|xXnhhhzyczy nn 则称欧拉方法稳定欧拉方法稳定。 注意注意: 这里 分别是以 为 初值得到的精确值，毫无舍入误 差，因此这里稳定性定义式对初 值的稳定性，即研究初值误差在 计算过程中的传递问题。 nn zy , 00, z y .21 定理定理1.4 1.4 在定理在定理1.21.2的条件下，欧拉方法是稳的条件下，欧拉方法是稳 定的。定的。 Th1.2 由定理1.2，我们看到如初始误差 ，则整体截断 误差的阶完全由局部截断误

11、差的阶决定，事实上，若局 部截断误差阶为 ，则整体截断误差阶为 。因 此为了提高数值算法的精度，往往从提高局部截断误差 的阶入手，这也时构造高精度差分方程数值方法的主要 依据。 0 0 )( 1p hO )( p hO .22 1.3 梯形法、隐式格式的迭代计算梯形法、隐式格式的迭代计算 dxxyxfdxyxf n n n n x x x x nn 11 )(,(),( )(,( nn xyxhf 在欧拉方法的推导过程，用矩形公式近似计算积分 1 ),()()( 1 n n x x nn dxyxfxyxy )(,( 2 )( 1 nn nn xyxf xx )(.( 11 nn xyxf 若

12、用梯形公式近似计算积分，则 .23 0 图图1.3 1n x n xx yxf, yxf, )(,( 2 )( ),( 1 1 nn nn x x xyxf xx dxyxf n n )(.( 11 nn xyxf .24 因此有 )(,()(,()( 2 1 )()( 1111 nnnnxnnn xyxfxyxfxxxyxy ),(),( 2 1 111 nnnnnn yxfyxfhyy(1.16) 这是一个隐式格式。 梯形公式局部截断误差： 首先给出后退欧拉公式的局部截断误差： )(,()()( 111nnnn xyxhfxyxy 后退欧拉公式 .25 假定 和解 充分光滑，则yxf, x

13、y * 111 )( nnn yxye )( 2 1 32 hOh 而欧拉公式的局部截断误差为 * 111 )( nnn yxye )( 2 1 32 hOh 梯形公式的局部截断误差为 )(,( 2 1 )()()( 1 1 * 11nnnn n nn xyxfhxyxyyxye ),( * 11 nn yxf )( 3 hO .26 类似于欧拉法，对梯形法和后退欧拉法也可平行 地建立它们的整体截断误差的阶分别为 和 ， 以及格式的收敛性和稳定性等定理。 )( 2 hO)(hO .27 ),(),( 2 111 nnnnnn yxfyxf h yy(1.16) 如何求解 ，采用迭代法，其格式如

14、下： 1n y 初始猜测 )0( 1 )( 11 )1( 1 ),(),( 2 n nn p nnn p n y yxfyxf h yy (1.18) 迭代法的收敛性 前已指出，梯形法是一个隐式格式 1 2 L h (1.19) 迭代法收敛的充分条件。 .28 当 ，有下面的预报-校正格式： 0p 也称为改进欧拉公式改进欧拉公式。 校正格式 预报格式 ),(),( 2 ),( )0( 111 )0( 1 nnnnnn nnnn yxfyxf h yy yxhfyy (1.20) .29 当然也可迭代多次： 校正格式 预报格式 ),(),( 2 ),( )( 11 )1( 1 )0( 1 p n

15、nnnn p n nnnn yxfyxf h yy yxhfyy (1.21) 当步长 取得适当小，用预报格式(欧拉法)已能算出 比较好的近似值，故迭代收敛很快，通常只需迭代二三次 就可满足精度要求，如果迭代多次仍不收敛，说明步长过 大，必须减少步长 ，再进行计算。 h h 梯形法较之欧拉法提高了精度，但增加了迭代次数， 因此增加了计算工作量。 .30 例例 1.21.2 试用预报校正格式(1.20)解初值问题 1| 1 , 01 0 x y xxyy (1.21) 1 . 0h取 。 从计算结果可以看出，欧拉法精度较低，预报 校正格式精度有所改善，大约精确到3为有效数 字。 .31 1.4

16、1.4 一般单步法、一般单步法、 格式格式KuttaRunge 前面，我们研究了欧拉法和梯形法，它们有 一个共同的特点，即在格式中只包括 的 值，或者说由 ，仅使用 的值计算出 的 值，这种格式称为单步格式，下面研究一般单步 法。 11, , nnnn yxyx 1 nn xx n y 1n y 1.4.1 1.4.1 一种构造单步法的方法一种构造单步法的方法泰勒级数法泰勒级数法 设初值问题 0 0 | ),( yy yxf dx dy xx .32 的解 阶可微，将 在 点展开为泰勒级 数，有 ) 1()(qxy xy 0 x )()( ! .)( ! 2 )()()( 1 0 )( 0 2

17、 000 qq q hOxy q h xy h xyhxyhxy (1.22) 由方程可得 ),()( 000 yxfxy yyxfyxfyxfxy yx ),(),(),()( 因此 ),(),(),()( 0000000 yxfyxfyxfxy yx )()(2)()( 2 yxyyyxyxx fffffffffxy dx d xy ,., xxyx fff 式中， 都是相对于变量的偏导数。于是式 .33 (1.22)可写成 )(),(,()()( 1 0000 q hOhxyxhxyhxy 其中 )(,( 1 1 1 1 00 00 | )(,( ! 1 (),(,( xyx q j j

18、 j j xyxf dx d h j hxyx 舍去 ，可得)( 1q hO ),( 0001 hyxhyy ),( 1112 hyxhyy ),( 1 hyxhyy nnnn 称 为一般单步法，显然局部截断误差 ),( 1 hyxhyy nnnn * )( nnn yxye .34 ),(,()()( 111 hxyxhxyxy nnnn q j n j j j nn xy dx d h j xyxy 1 11 )( ! 1 )()( )( 1 q hO 所以局部截断误差为 ，在式(1.23)中令 ， 即得欧拉法。 )( 1q hO 1q .35 为任意关于 的函数，其对于 微分方程 ),(

19、,(hxyx),(,(hxyx ),(yxf dx dy 的解 满足 )(xy )(),(,()()( 1 q hOhxyxhxyhxy (1.24) 且 为使上式成立的最大整数，则称 q ,.2 , 1 , 0),( 1 nhyxhyy nnnn (1.25) 为 阶单步法，欧拉法为一阶单步法，泰勒级数 法式(1.23)为 阶单步法。 q q 定义定义 1.2 1.2 给出单步法 ),( 1 hyxhyy nnnn 1.4.2 1.4.2 一般单步法基本理论一般单步法基本理论 .36 定义定义 1.31.3 如果， 则称单步法 为与初值问题(1.3)相容的。 ),()0 ,(yxfyx ),

20、( 1 hyxhyy nnnn )(),()()( 1 q hOhyxhxyhxy 从而 )0 ,( )()( lim 0 yx h xyhxy h 由单步法的定义得: 因此有如下定义: .37 定理定理 1.51.5 如果 对于 ， 以及所有实数 满足 条件，则单步 法(1.23)稳定。 满足：的解 nn zy , ),(hyx Xxx 0 0 0hh y Lipschitz | 00 zyczy nn 0 1 ),( z hzxhzz nnnn 0 1 ),( y hyxhyy nnnn 欲使定理成立要证明存在常数 , 对 0, 0 0 hc .38 定理定理 1.61.6 如果 对于 ，

21、 以及所有实数 关于 满足 条件，则 收敛的充要条件是格式相容，即 满足 。 ),(hyx Xxx 00 0hh y hyx, Lipschitz ),( 1 hyxhyy nnnn ),()0 ,(yxfyx 由相容性可以得到格式的收敛性: .39 定理定理 1.71.7 在定理1.5的条件下，如果局部截断 误差 为 ，则单步法 的整体截断误差 满足 n n R)( 1q hO ),( 1 hyxhyy nnnn nnn yxy)( ) 1(| )( 0 )( 001 xXLqxXL n e L c heee(1.26) 特别若 ,则 ，整体截断误差比局部截 断误差低一阶。 0 0 )( q

22、 n hO 关于单步法的整体截断误差 ，有: .40 1.4.3 1.4.3 格式格式KuttaRunge 从前面讨论可见，构造高阶单步法的关键在 于构造 ，使),(hyx )(),(,()()( 1 1 q nnnn hOhxyxhxyxy 中的局部截断误差阶尽可能高。前面我们利用泰 勒级数法构造了一个欧拉法，这时 ，局部截断误差与 同阶，这是一个一阶格式。 为了要求 ，利用泰勒级数法得到一个二阶格 1),(),(,()()( 2 1 qhOhxyxhfxyxy nnnn )( 2 hO 2q ),(),(),( 2 ),( 2 1nnnnynnxnnnn yxfyxfyxf h yxhfy

23、y (1.27) .41 这时我们有 )(,()()( 1nnnn xyxhfxyxy )()(,()(,()(,( 2 3 2 hOxyxfxyxfxyxf h nnnnynnx 格式(1.27)计算过程中要求函数 的二个 偏导数 在 处的值，比较麻烦，可以预 计，利用泰勒级数法推导出的高阶格式需要求更 多的偏导数值，计算繁复。那么是否可以避免计 算偏导数，而得到高阶单步格式的 呢？分 析梯形法的预报校正格式(1.20) ),(yxf yx ff ,)(,( nn xyx ),(hyx 校正格式 预报公式 ),(),( 2 ),( )0( 111 )0( 1 nnnnnn nnnn yxfy

24、xf h yy yxhfyy .42 二级二阶 方法。一般而言，二级二 阶 格式可以写成KuttaRunge KuttaRunge ),( ),( )( 12122 1 22111 KhbyhaxfK yxfK KcKchyy nn nn nn (1.28) 适当选择参数 ，使局部截断误 差 21221 ,bacc )(,()()( 111nnnnn xyxfchxyxyR )(,()(,( 2122nxnn xyxfhbxyhaxfc )( 3 hO .43 由 )(,( 2 )(,()( 2 1nnxnnnn xyxf h xyxhfxyR )() )(,()(,( 3 hOxyxfxyx

25、f nnnny )(,()(,()( 21nnnnn xyxfhcxyxfhcxy )() )(,()(,()(,( 3 212 hOxyxfhbxyxfhaxyxf nnnnynnx 因此要求满足 2 1 2 1 1 221 22 21 cb ca cc (1.29) .44 (1) 取 ，则 ,即得二级二阶 法 2 1 1 c 1, 2 1 2122 bac KuttaRunge ),( ),( )( 2 12 1 211 hKyhxfK yxfK KK h yy nn nn nn (1.30) 0 1 c 2 1 , 1 2122 bac ) 2 1 , 2 1 ( ),( 12 1 2

26、1 hKyhxfK yxfK hKyy nn nn nn (1.31) 这是一个含有四个参数、三个方程的方程组,因此 由一个自由参数，解答不唯一。 改进欧拉公式 变形欧拉公式 (2)令 ， ，由此得算式为 .45 (3) 取 ， ，则有 4 1 1 c 3 2 , 4 3 2122 bac ) 3 2 , 3 2 ( ),( )3( 4 12 1 211 hKyhxfK yxfK KK h yy nn nn nn (1.32) 根据同样的思想,可以构造更高阶精度的 Runge-Kutta 方法。 ),( ),( ),( )( 23213133 12122 1 3322111 hKbhKbyha

27、xfK hKbyhaxfK yxfK KcKcKchyy nn nn nn nn (1.33) 三级三阶 Runge-Kutta法一般可以写成： .46 6 1 6 1 3 1 )( 3 1 )( 3 1 2 1 )( 2 1 1 21323 3223 2 32313 2 212 3231332122 2 33 2 22 32313212 3322 32 bbc bac bbcbc bbacbac acac bbcbc acac ccc (1.34) )( 4 1 hORn 故要求 ，必须有 .47 (1) 令 ，则 ， 6 1 31 cc1, 2 1 , 6 4 322 aac1, 2, 2

28、 1 313221 bbb 故有 三级三阶算法Kutta )2,( ) 2 , 2 ( ),( )4( 6 213 1 2 1 3211 hKhKyhxfK hK y h xfK yxfK KKK h yy nn nn nn nn (1.36) 令 ，解得 ， 3 2 , 3 1 32 aa 4 3 , 0, 3 1 3221 ccb 0, 3 2 , 4 1 31321 bbc ，故有 三级三阶 算法HeumKR .48 ) 3 2 , 3 2 ( ) 3 , 3 ( ),( )3( 4 213 1 2 1 311 hKhKyhxfK hK y h xfK yxfK KK h yy nn n

29、n nn nn (1.37) 式(1.36)，(1.37)是三级三阶 格式。 KuttaRunge 格式的局部截断误差为 。)( 4 hO .49 可以设计四级四阶 格式KuttaRunge ),( ),( ),( ),( )( 34324214144 23213133 12122 1 443322111 hKbhKbhKbyhaxfK hKbhKbyhaxfK hKbyhaxfK yxfK KcKcKcKchyy nn nn nn nn nn (1.38) .50 为了达到四级四阶格式，可得13个参数满足 11个方程 24 1 12 1 )( 8 1 )( 6 1 )( 4 1 3 1 2

30、1 1 433224 43 2 342 2 2432 2 23 4433422432323 43342243223 3 44 3 33 3 22 2 44 2 33 2 22 443322 4321 4342414 32313 212 bbac babacbac ababacbaac babacbac acacac acacac acacac cccc bbba bba ba .51 (1) 经典四级四阶 格式取定 ， 则得： KuttaRunge 2 1 32 aa ),( ) 2 , 2 ( ) 2 , 2 ( ),( )122( 6 34 2 3 1 2 1 43211 hKyhxfK

31、hK y h xfK hK y h xfK yxfK KKKK h yy nn nn nn nn nn (1.39) 这是最为著名的经典四级四阶 格式。KuttaRunge 格式的局部截断误差为 。) ( 5 hO .52 1.5 1.5 线性多步法线性多步法 前面利用欧拉法或梯形法求未知函数在 的近似值 ，基本思想是在积分表达式 1 n xx 1n y dxxyxfxyxy n n x x nn 1 )(,()()( 1 中被积函数 ，用水平直线 或连接 两点的直线 代替。然而为了近 似 中的曲线 也可用多点插值曲线， 如我们利用 插值可得经过 的曲线 为 1 ,),( nn xxxyxf)

32、(,(),( nn xyxfyxf )(,(,( nnn xyxfx)(,(,( 111nnn xyxfx )(,(),( 1 nn n xyxf h xx yxf)(,( 11 nn n xyxf h xx 1 , nn xx )(,(xyxf Lagrange )(,(,( 222nnn xyxfx )(,(,( , 111nnn xyxfx )(,( , nn xyx )(,( 2 xyxL )(,( 2 )( )(,( 22 2 1 2nn nn xyxf h xxxx xyxL .53 )(,( )( 11 2 2 nn nn xyxf h xxxx )(,( 2 )( 2 12 n

33、n nn xyxf h xxxx (1.50) 用它近似 中的 ，则得到积分近似值 1 , nn xx)(,(xyxf dxxyxLxyxy n n x x nn 1 )(,()()( 21 )(,( 12 5 )( 22nnn xyxfhxy )(,( 12 23 )(,( 12 16 11nnnn xyxfxyxf 像欧拉格式一样，我们可得近似求解格式如下： ),( 12 5 221nnnn yxfhyy ),( 12 23 ),( 12 16 11nnnn yxfyxf (1.51) .54 与欧拉格式 不同之处在于增加了包 括 , 的两项,即由 计算 ，与欧拉格式仅由前面一点 计算 的

34、这种单步法不同，格式(1.51)称为多步法。 现在，设已给出常微分方程初值问题 ),( 1nnnn yxhfyy ),( 22nn yxf ),( 11nn yxf ),(),(),( 1122nnnnnn yxyxyx 1n y ),( nn yx 1n y 0 0 | ),( yy yxf dx dy xx 的解 在 处的近似值 ，或者说 给出表头 )(xy n xxx,., 10 n yyy,., 10 x y n xxxx 210 n yyyy 210 .55 研究如何由表头给出 处的近似 值 。 根据 Xhnxxn ) 1( 01 1n y 1 )()()( 1 n n x x nn

35、 dxxyxyxy (1.52) 和 ，利用表头的值和 插值法求 出 的近似表达式，再利用(1.52)就得到 的 近似值。 具体操作如下： 用 表示用 处 的 值构造出的 的 插值多项式，用 表示 相应的插值余项，即 从而 ),()(yxfxyLagrange )(x y )( 1n xy )( , xL kn knnn xxx ,., 1 )(),.,(),( 1knnn xyxyxy )(x y Lagrange kn r , )()()( , xrxLxy knkn 11 )()()()( ,1 n n n n x x kn x x knnn dxxrdxxLxyxy .56 舍去余项

36、，并用 代替 ，则可得 的近似值的表达式 1 )( , n n x x knkn dxxrR j y)( j xy)( 1n xy dxLyy n n x x knnn 1 * ,1 (1.53) 为 中 的值用 代替， 为局部截断误差。 上述格式中，被插值点 不包括在插值 基点所决定的最大区间 内，故称为外插公式。 由常微分方程， 根据 ， ， ， 得插值公式 为 ，由于插值基点为等距，令 ，利用牛顿后插公式，有 )( * , xL kn )( , xL kn )( j xy j ykn R , )( 1 nn xxxx nkn xx, )(,()(xyxfxy nnn yxfx, 111 , nnn yxfx knknkn yxfx , )( * , xL kn 10 ,tthxx n . ! 2 ) 1( ! 1 )( 2 2 1