十二节常微分方程数值解法

1、12 欧拉（欧拉（euler）方法、向后欧拉法、梯形法及梯形）方法、向后欧拉法、梯形法及梯形法的预估校正法法的预估校正法 欧拉法的收敛性欧拉法的收敛性 龙格库塔方法、线性多步法、预估校正法龙格库塔方法、线性多步法、预估校正法*。 一阶微分方程组与高阶微分方程的数值解法一阶微分方程组与高阶微分方程的数值解法*3在解决科技领域的实际应用问题时，常微分方程求解在解决科技领域的实际应用问题时，常微分方程求解是常见的。本章着重讨论一阶方程初值问题是常见的。本章着重讨论一阶方程初值问题的数值解法。对高阶方程和微分方程组的数值解，的数值解法。对高阶方程和微分方程组的数值解，其基本思想是完全一样的其基本思想是

2、完全一样的解初值问题有多种解解初值问题有多种解析方法，但解析法只能对一些特殊类型的方程才析方法，但解析法只能对一些特殊类型的方程才能求出其准确解，多数情况只能用近似方法求解。能求出其准确解，多数情况只能用近似方法求解。初值问题的数值解法，就是寻求方程的解初值问题的数值解法，就是寻求方程的解( )y x在自变量在自变量x的一系列离散节点上的近似值。的一系列离散节点上的近似值。00(,),(),dyfxydxy xy4,)()(),(3212100/nnyyyyxxxxyyxyyxfy近似解：，处的在节点求：精确解初值问题初值问题5 相邻两节点间的距离相邻两节点间的距离 称为步长，通称为步长，通常

3、在计算上采用相等的步长常在计算上采用相等的步长 ，这时等距，这时等距节点节点 ， 初值问题的数值解法的基本特点是：求解过程是初值问题的数值解法的基本特点是：求解过程是顺着节点排列的顺序一步一步的向前推进，即按顺着节点排列的顺序一步一步的向前推进，即按递推方法由已知的递推方法由已知的 求出求出 。所。所以，初值问题的数值解法就是建立这种递推公式。以，初值问题的数值解法就是建立这种递推公式。1nnnhxxnhh0nxxnh0,1,2,n 01,nyyy1ny6将微分方程两端从将微分方程两端从nx到到1nx积分，得积分，得11()( )( , ( )nnxnnxy xy xf x y x dx(0,

4、1,2,)n这样，求原初值问题式的解，转化为求问题式这样，求原初值问题式的解，转化为求问题式的解的解, ,利用各种求积公式就可以得到一些求利用各种求积公式就可以得到一些求()ny x的近似公式。的近似公式。 7 差商方法差商方法001001/)(),()()()(),(yxyyxhfyyyxyhxyxyyxfynnnnnn8 数值积分方法数值积分方法)(,()(,()(,()()()(,()()(,)(,(/nnhxxhxxnnnhxxxyxhfdttytfdttytfxyhxyxxdttytfxyhxyhxxxyxfynnnn时，有当得上积分，在9 数值积分方法数值积分方法100()()(

5、, ( )(, ()(,)()nnxhnnxnnnnnny xhy xf t y t dthf xy xyyhf xyyy x（看成矩形）10 向后差商向后差商00111001/)(),()()()(),(yxyyxhfyyyxyhxyxyyxfynnnnnn11 中心差商中心差商00110011/)(),()(2)()(),(yxyyxhfyyyxyhxyxyyxfynnnnnn12)(),(),(2)(,()(,(2)(,()()(0011111xyyyxfyxfhyyxyxfxyxfhdttytfxyhxynnnnnnnnnnhxxnnnn13 梯形公式梯形公式( (见上页见上页) )，

6、实际上是，实际上是euler方法和隐式方法和隐式euler方法的算术平均。方法的算术平均。 梯形公式的精度为二阶。梯形公式的精度为二阶。 例：用梯形公式求下列初值问题的解在例：用梯形公式求下列初值问题的解在 ).01. 0(01. 0yx上的值1)0(,yydxdy14 改进的改进的euler方法为方法为euler方法和梯形公式的结合，方法和梯形公式的结合，也称作预估也称作预估-校正法。校正法。)(),(),(2),(001111xyyyxfyxfhyyyxhfyynnnnnnnnnn15 嵌套形式嵌套形式),(,(),(211nnnnnnnnyxhfyxfyxfhyy1ny16),(),(,

7、(),(),()(21111121211hkyxfyxhfyxfyxfkyxfkkkhyynnnnnnnnnnnn平均化形式17 称一种数值方法是称一种数值方法是p p阶的，如果其局部截断误差阶的，如果其局部截断误差为为 。 eulereuler方法和隐式方法和隐式eulereuler方法的精度是一阶的。方法的精度是一阶的。 二步二步eulereuler方法的精度是二阶的。方法的精度是二阶的。)(1pho18 改进的改进的eulereuler方法也可写成方法也可写成),(),(22121211kyhxhfkyxhfkkkyynnnnnn19),(),(12122122111kbyhaxhfky

8、xhfkkckcyynnnnnn20 要使二阶方法的局部截断误差为要使二阶方法的局部截断误差为 ，四个系，四个系数值应满足下列关系式：数值应满足下列关系式：1;21;12212221abaccc)(3ho21 特例特例1 1： 2122111212111;121122(,)(,)nnnnnnaccbyykkkhfxykhfxhyk令， 则这 就 是 前 面 讲 的 改 进 的 尤 拉 法 。22 特例特例2 2：称为变形的尤拉法。该式称为中点方法，也，令)21,21(),(21;211,01212121221kyhxhfkyxhfkkyybaccnnnnnn231123121312 3 (4)

9、6(,)(,)22(,2)nnnnnnnnprkrkhyykkkkf xyhhkf xykkf xh yhkhk类似地，对，即三个点，通过更复杂的计算，可导出三阶公式。常用的三阶公式为：三阶龙格库塔方法24四阶龙格库塔方法112341213243 4 (22)6(,)(,)22(,)22(,)nnnnnnnnnnprkrkhyykkkkkf xyhhkf xykhhkf xykkf xh yhk对，即四个点，可导出四阶公式。常用的四阶公式为： 25112341211322433 h = 0 .2 ,x = 0 x = 12( 01);( 0 )1 .(22);62;2;222;222 ().n

10、nnnnnnnnnnnnnxyyxyyhyykkkkxkyyxhhkykhykxhhkykhykxhkyh kyh k例 ： 设 取 步 长从直 到用 四 阶 龙 格 库 塔方 法 求 解 初 值 问 题解 ： 由 经 典 的 四 阶 龙 格 库 塔 公 式 得261rkrk46rk52rk rkrk rk两 点 说 明 ：） 当 p=1,2,3,4时 ，公 式 的 最 高 阶 数 恰 好 是 p,当 p4时 ，公 式 的 最 高 阶 数 不 是 p， 如 p=5时 仍为 ， p= 时公 式 的 最 高 阶 数 为 。）方 法 的 导 出 基 于 taylor展 开 ， 故 要 求 所 求 问

11、题 的 解 具 有 较 高 的 光 滑 度 。当 解 充 分 光 滑 时 ， 四 阶方 法 确 实 优 于 改 进euler法 。 对 一 般 实 际 问 题 ， 四 阶方 法 一 般 可 达到 精 度 要 求 。如 果 解 的 光 滑 性 差 ， 则 用 四 阶方 法 解 的 效 果不 如 改 进 euler法 。27()1()51115(2 )15(2 )11(2 )11()11,(),2,2()2,2()1.()1 6hnhnnnnhnhnnhnnhnnhyyxyc hhxxhychyxycyxyyxyn以 经 典 四 阶 龙 格 库 塔 公 式 为 例 。 从 节 点 x 出 发 ，

12、以为步 长 求 一 近 似 值将 步 长 折 半 ， 即 取为 步 长 从跨 两 步 到， 求 一 近 似值每 跨 一 步 的 截 断 误 差 是因 此 有由 上 两 式 (2 )(2 )()11111().1 5hhhnnnnyxyyy变步长的龙格库塔方法2811111r k(,)(,)(2 , 3,)(,)()pnniiinniininijjjnnnyyhc kkfxykfxa hyhb kipxyt a y lo ryxxt a y lo r一 般 地 ，方 法 设 近 似 公 式 为确 定 原 则 是 使 近 似 公 式 在处 的展 开 式 与在处 的展 开 式 的 前 面 项 尽 可

13、 能 多 地 重 合 。111112121(,)11()22(,)(,)nnnnnnnnnnyyh ke u le rkfxyyyhkke u le rkfxykfxhyh k公 式 改 进公 式2911-r1rk,nnnnnnyyyyyy单步法在计算时，只用到前一步的信息 。为提高精度，需重新计算多个点处的函数值，如方法，计算量较大。如何通过较多地利用前面的已知信息，如 ，来构造高精度的算法计算，这就是多步法的基本思想。11110111 (,),(,) ,(,) (1, ,)00 taylonnn rnnn rn rrrnin iin iiiiikkkyyyf xyf xyyyhfff xy

14、knnnr多步法中最常用的是线性多步法，它的计算公式中只出现， ，及的一次项，其一般形式为其中均为常数，。若 ，显式；，隐式。构造线性多步公式常用r展开和数值积分方法。线性多步法30nn taylorxtaylor)xtaylor,iin+1利 用展 开 导 出 的 基 本 方 法 是 ： 将 线 性 多 步 公 式在处 进 行展 开 ， 然 后 与 y(x在处 的展开 式 相 比 较 ， 要 求 它 们 前 面 的 项 重 合 ， 由 此 确 定 参 数。101111011 1( ) ()nnnnnnry xyyyhfff以为例：设初值问题的解充分光滑，待定的两步公式为()()()21()

15、(1, 2,),( )( )()()()2!()kknnnppnnnnnnnpnyyxky xxtayloryyy xyyxxxxxxpoxx记则在处 的展 开 为线性多步公式的导出31231( 4 )( 5 )45( 6 )1111() ,()(,) () , ()2 !3 ! ()4 !5 ! (,)()iiiiinnnnnnnnnnnnnnnyyxyxfxyinyyyyxhyyhhhyyhhohffxyyxyyh假 设 前步 计 算 结 果 都 是 准 确 的 ， 即则 有2( 4 )( 5 )34( 5 )21111( 4 )(32 ! ()3 !4 ! (,) (,)()2 ! 3

16、!nnnnnnnnnnnnnnnnyhyyhhohffxyyyffxyyxyyhhyyh5 )4( 5 )()4 !hoh线性多步公式的导出32211011101113(4)4111111(5)56111()()()2()()6222466()()1202424nnnnnnnyyy hy hy hyhy ho h 将以上各公式代入并整理，得1(5 )2561()()() 2!5!p+1 nnnnnnnpy xxtayloryyy xyy hhho h为 使 上 式 有阶 精 度 ， 只 须 使 其 与在处 的展 开 式的 前项 重 合 。线性多步公式的导出3301010111111111111

17、11221111622611112 4662 4aaaaaa5,5,1iip个参数 只须 个条件。由推导知，如果选取参数，使其满足前个方程（p=1,2,3,4)，则近似公式为p阶公式。线性多步公式的导出3411011111,0,02 ()2nnnnhyyff0如满 足 方 程 组 前 三 个方 程 ， 故 公 式此 为 二 阶 公 式 。0111011115(5)61140,1,33 (4)31 ()90nnnnnnnhyyfffrh yo h 又如：解上面方程组得相应的线性二步四阶公式（simpson公式)为其截断误差为由此可知，线性二步公式至多是四阶公式。线性多步公式的导出35101011

18、11 23()1()()1(1,2,) 1()(1)!nnnriirrkkiiiippnirxtaylory xxtaylorikikpphrip一般地，线性多步公式中有个待定参数，如令其右端在处的展开式与在处的展开式的前p+1项系数对应相等，可得方程组其解所对应的公式具有 阶精度，局部截断误差为(1)1121)() ()rrppiniippiyo h显然，线性多步公式至多可达到2r+2阶精度。线性多步公式3612310101001231123 (adams)r=30,1()()1(1, 2,3 4)5559379=1,24242424(555937924riirrkkiiiinnnnnnik

19、ikhyyffff （ 一 ） 阿 达 姆 斯公 式取， 并 令由 方 程 组，可 解 得，相 应 的 线 性 多 步 公 式 为1)=0adams因， 此 式 称 为显 式 公 式 ， 是 四 阶 公 式 .常用的线性多步公式3753354(5 )61115(5 )61()5()()5!251()720niiniinhriiyo hh yo h 局 部 截 断 误 差 为12330101211125( 5 )610,91 951 = 1,2 42 42 42 4(91 95)2 4a d am s1 9()7 2 0nnnnnnnnhyyffffrhyoh 如 果 令由 方 程 组 可 解

20、得，相 应 的 线 性 多 步 公 式 为称 其 为 四 阶隐 式 公 式 ， 其 局 部 截 断 误 差 为常用的线性多步公式381111+11()()( ,( )( ),( )( )( ),1nnnnxxnnxxnnnrnnnry xy xf x y xdxf x dxxxxxxxf xxf xr r基本思想是首先将初值问题化成等价的积分形式用过节点或的的r次插值多项式代替求积分 即得阶的线性多步公式。123330123303,( ) ( )( )()()()()()( )(0,1,2,3)()()nnnniniinnnninininjjjirxxxxf xlxlx f xxxxxxxxx

21、lxixxxx例如时，过节点的三次插值多项式为其中利用数值积分方法求线性多步公式391111131301233231313233()()( )( )()()()()()6()()()()2()()()()2()()nnnnnnnnnnxxnninixxixnnnnxxnnnnxxnnnnxnny xy xlx dxlx dx f xxxxxxxf xdxhxxxxxxf xdxhxxxxxxf xdxhxxxf x1123123)()655()59()37()9()24nnxnnxnnnnxxxdxhhf xf xf xf x利用数值积分方法求线性多步公式401111233,(),(),(,)

22、()(,() (,1,2,3),(5559379)24,nnnnkkkkkknnnnnnnnyyy xy xfxyfxfxy xkn nnnhyyffffadamsxxadams对 上 式 用代 替用代 替则 得这 就 是 四 阶显 式 公 式 。 由 于 积 分 区 间 在 插 值区 间外 面 ， 又 称 为 四 阶外 插 公 式 。111(4)(5)3310031(5)35(5)10 ()()()()4!4!,),( )251()( )4!720nnnnnnxxxxnnjnjxxjjnnxnnjxjfyrxxdxxxdxxxyrxxdxh y由插值余项公式可得其局部截断误差为由积分中值定理，存在（使得利用数值积分方法求线性多步公式41112231112311 ,( ) ( )( )()()()()()( )(1,0,1, 2)()()( ) nnnniniinnnninininjjjinnxxxxf xlxlx