1、1第第9 9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法 本章着重考察一阶方程的初值问题.)(),(00yxyyxfy（1.1）（1.2）只要函数 适当光滑譬如关于 满足利普希茨(Lipschitz）条件),(yxfy.),(),(yyLyxfyxf（1.3）理论上就可以保证初值问题（1.1），（1.2）的解 存在并且唯一.)(xyy 2 所谓数值解法数值解法，就是寻求解 在一系列离散节点)(xy 121nnxxxx上的近似值 .相邻两个节点的间距 称为步长步长.,121nnyyyynnnxxh 1 如不特别说明，总是假定 为定数，这时节点为 .),2,1(ihhi),2,1,0
2、(0inhxxn 首先，要对方程（1.1）离散化，建立求数值解的递推公式.一类是计算 时只用到前一点的值 ，称为单步单步法法.另一类是用到 前面 点的值 ，称为 步法步法.1nyny1nyk11,knnnyyyk 其次，要研究公式的局部截断误差和阶，数值解 与精确解 的误差估计及收敛性，还有递推公式的计算稳定性等问题.ny)(nxy3 9.1 9.1 简单的数值方法简单的数值方法.)(),(00yxyyxfy一 欧拉法与后退的欧拉法（1.1）（1.2）如果对方程（1.1）从 到 积分，得 nx1nx.)(,()()(11nnxxnndxxyxfxyxy（1.4）右端积分用左矩形公式 近似，再以
3、 代替 代替 得到欧拉公式)(,(nnxyxfhny1),(nnyxy)(1nxy).,(1nnnnyxhfyy（1.5）4图8-1若初值 已知，则依公式（1.5）可逐步算出),(),(11120001yxhfyyyxhfyy0y欧拉法又称折线法，几何意义如下图 5 解解 欧拉公式的具体形式为).2(1nnnnnyxyhyy取步长 ，计算结果见表9-1.1.0h7321.17848.10.14142.14351.15.06733.17178.19.03416.13582.14.06125.16498.18.02649.12774.13.05492.15803.17.01832.11918.12
4、.04832.15090.16.00954.11000.11.0)()(19nnnnnnxyyxxyyx计算结果对比表 例例1 1 求解初值问题.1)0(),10(2yxyxyyxy21精确解：6 如果在（1.4）中右端积分用右矩形公式 近似，则得另一个公式)(,(11nnxyxfh),(111nnnnyxhfyy（1.6）称为后退的欧拉法后退的欧拉法.后退的欧拉公式与欧拉公式有着本质的区别，后者是关于 的一个直接的计算公式，这类公式称作是显式的显式的；1ny（1.4）然而公式（1.6）的右端含有未知的 ，它实际上是关于 的一个函数方程，这类公式称作是隐式的隐式的.1ny1ny1)(,()()
5、(1nnxxnndxxyxfxyxy7 隐式方程（1.6）通常用迭代法求解.先用欧拉公式),()0(1nnnnyxhfyy给出迭代初值 ，用它代入（1.6）式的右端，使之转化为显式，直接计算得)0(1ny),()0(11)1(1nnnnyxhfyy然后再用 代入（1.5）式，又有)1(1ny).,()1(11)2(1nnnnyxhfyy8如此反复进行，得).,1,0(),()(11)1(1kyxhfyyknnnkn（1.7）由于 对 满足利普希茨条件（1.3）.由（1.7）减（1.6）得),(yxfy.),(),(1)(111)(111)1(1nknnnknnnknyyhLyxfyxfhyy由
6、此可知，只要 迭代法（1.7）就收敛到解 .1hL1ny9二二 梯形方法梯形方法 为得到比欧拉法精度高的计算公式，在等式（1.4）右端积分中若用梯形求积公式近似，并用 代替 代替 ，则得 ny1),(nnyxy)(1nxy),(),(2111nnnnnnyxfyxfhyy（1.8）称为梯形方法梯形方法.梯形方法是隐式单步法，可用迭代法求解.同后退的欧拉方法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为 1)(,()()(1nnxxnndxxyxfxyxy（1.4）10).,2,1,0(),(),(2);,()(11)1(1)0(1kyxfyxfhyyyxfhyyknnnnnknnnnn（
7、1.9）为了分析迭代过程的收敛性，将（1.8）式与（1.9）式相减，得(1)()111111(,)(,),2kknnnnnnhyyfxyfxy于是有,2)(11)1(11knnknnyyhLyy式中 为 关于 的利普希茨常数.),(yxfyL11如果选取 充分小，使得 h,12hL则当 时有 ，这说明迭代过程（1.9）是收敛的.k1)(1nknyy12三三 单步法的局部截断误差与阶单步法的局部截断误差与阶 初值问题(1.1)，(1.2)的单步法可用一般形式表示为),(11hyyxhyynnnnn（1.10）其中多元函数 与 有关，当 含有 时，方法是隐式的，若不含 则为显式方法，所以显式单步法
8、可表示为),(yxf1ny1ny),(1hyxhyynnnn（1.11）称为增量函数，例如对欧拉法（1.5）有),(hyx).,(),(yxfhyx13 对一般显式单步法则可如下定义.定义定义1 1 设 是初值问题(1.1)，(1.2)的准确解，称)(xyy),(,()()(11hxyxhxyxyTnnnnn（1.12）为显式单步法（1.11）的局部截断误差局部截断误差.之所以称为局部的，是假设在 前各步没有误差.1nTnx当 时，计算一步，则有)(nnxyy.),(,()()(),()()(11111nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy所以，局部截断误差可理解为用方
9、法（1.11）计算一步的误差.14 根据定义，显然欧拉法的局部截断误差),()(2)()()()(,()()(3211hOxyhxyhxyhxyxyxhfxyxyTnnnnnnnnn 这里 称为局部截断误差主项.显然 .)(22nxyh)(21hOTn15 定义定义2 2 设 是初值问题(1.1)，(1.2)的准确解，若存在最大整数 使显式单步法(1.11)的局部截断误差满足)(xyp),(),()()(11pnhOhyxhxyhxyT（1.13）则称方法（1.11）具有 阶精度阶精度.p 以上定义对隐式单步法（1.10）也是适用的.16 对后退欧拉法（1.6）其局部截断误差为).()(2)(
10、)()()()(2)()(,()()(322321111hOxyhhOxyhxyhhOxyhxyhxyxhfxyxyTnnnnnnnnnn 这里 ，是1阶方法，局部截断误差主项为 .1p)(22nxyh 17 对梯形法（1.8）有).()(12)()(2)()()(2)(!3)(2)()()(2)()(434232111hOxyhhOxyhxyhxyxyhxyhxyhxyhxyxyhxyxyTnnnnnnnnnnnnn 所以梯形方法（2.7）是二阶的，其局部误差主项为 .)(123nxyh 18四四 改进的欧拉公式改进的欧拉公式 梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式（1.9）进
11、行实际计算时，每迭代一次，都要重新计算函数 的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测.),(yxf 为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 ,称之为预测值预测值，预测值 的精度可能很差，再用梯形公式（1.8）将它校正一次，即按（1.9）式迭代一次得，这个结果称校正值校正值，而这样建立的预测-校正系统通常称为改进改进的欧拉公式：的欧拉公式：1ny1ny19预测),(1nnnnyxhfyy校正).,(),(2111nnnnnnyxfyxfhyy（1.14）或表示为下列平均化形式).(21),(),(11cpn
12、pnncnnnpyyyyxhfyyyxhfyy20 例例2 2 用改进的欧拉方法求解初值问题 解解 改进的欧拉公式为).(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy仍取 ，计算结果见表9-2.同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.1.0h.1)0(),10(2yxyxyy217321.17379.10.14142.14164.15.06733.16782.19.03416.13434.14.06125.16153.18.02649.12662.13.05492.15525.17.01832.11841.12.04832.14860.16.
13、00954.10959.11.0)()(29nnnnnnxyyxxyyx欧拉法）计算结果对比（改进的表7321.17848.10.14142.14351.15.06733.17178.19.03416.13582.14.06125.16498.18.02649.12774.13.05492.15803.17.01832.11918.12.04832.15090.16.00954.11000.11.0)()(19nnnnnnxyyxxyyx）计算结果对比（欧拉法表22 9.3 9.3 龙格龙格-库塔方法库塔方法 ).,(),(2111nnnnnnyxfyxfhyy9.3.1 9.3.1 显式龙格
14、显式龙格-库塔法的一般形式库塔法的一般形式 上节给出了显式单步法的表达式),(1hyxhyynnnn其局部截断误差为(前提：)11111()()()(,(),)()pnnnnnnnTy xyy xy xhxy xhO h对欧拉法 ，即方法为 阶，若用改进欧拉法，它可表为)(21hOTn1p).,(,(),(21nnnnnnnnyxhfyhxfyxfhyy（3.1）),(1nnnnyxhfyy()nnyy x23此时增量函数).,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx（3.2）它比欧拉法的 ，增加了计算一个右函数 的值，可望 .),(),(nnnnyxfhyxf2p
15、若要使得到的公式阶数 更大，就必须包含更多的 值.pf 从方程（1.1）等价的积分形式（1.4），即,)(,()()(11nnxxnndxxyxfxyxy（3.3）若要使公式阶数提高，就必须使右端积分的数值求积公式24精度提高，必然要增加求积节点，为此可将（3.3）的右端用求积公式表示为1.)(,()(,(1nnxxriininihxyhxfchdxxyxf点数 越多，精度越高，上式右端相当于增量函数 ，为得到便于计算的显式方法，可类似于改进欧拉法（3.1），（3.2），将公式表示为 r),(hyx),(1hyxhyynnnn（3.4）其中,2),(),(,),(1111riKhyhxfKyx
16、fKKchyxijjijnininnriiinn（3.5）25这里 均为常数.ijiic,(3.4)和(3.5)称为 级显式龙格龙格-库塔库塔（Runge-Kutta）法法，简称R-K方法.r 当 时，就是欧拉法，此时方法的阶为 .),(),(,1nnnnyxfhyxr1p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，下面将证明阶 .2r2p 要使公式(3.4),(3.5)具有更高的阶 ,就要增加点数 .pr 下面只就 推导R-K方法.2r26 9.3.2 9.3.2 二阶显式二阶显式R-KR-K方法方法 对 的R-K方法，由（3.4），（3.5）可得到如下的计算公式 2r).,(
17、),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6）这里 均为待定常数，希望适当选取这些系数，使公式阶数 尽量高.21221,ccp 根据局部截断误差定义，（3.6）的局部截断误差为 111112221()()()(,)(,),nnnnnnnnnnTy xyy xy xh c fxycfxh yhf（3.7）27这里 .),(),(nnnnnyxffxyy 为得到 的阶 ，要将上式各项在 处做泰勒展开，由于 是二元函数，故要用到二元泰勒展开，各项展开式为 1nTp),(nnyx),(yxf),(!32)(4321hOyhyhyhyxynnnnn 其中 )
18、;,(),(),(),(),(2),(,),(),()(,(,),(2nnynnnxnnynnyynnnxynnnxxnnnnynnxnnnnnnnyxffyxfyxfyxffyxffyxfyfyxfyxfxyxfdxdyfyxfy（3.8）28).(),(),(),(2212212hOfhyxfhyxfffhyhxfnnnynnxnnnn将以上结果代入（3.7）则有).(),()21(),()21()1()(),(),(),(),(2322122222132122121hOhfyxfchyxfchfcchOhfyxfhyxffcfchfyxfyxfhfhTnnnynnxnnnnynnxnnn
19、nnynnxnn要使公式（3.6）具有 阶，必须使 2p29.021,021,012122221cccc（3.9）即.1,21,212121222cccc（3.9）的解是不唯一的.令 ，则得 02 ac.21,12121aac这样得到的公式称为二阶R-K方法，如取 ，则2/1a.1,2/121221cc这就是改进欧拉法（3.1）.30若取 ，则 .得计算公式 1a2/1,1,021221cc).2,2(),(,12121KhyhxfKyxfKhKyynnnnnn（3.10）称为中点公式中点公式，相当于数值积分的中矩形公式.（3.10）也可表示为).,(2,2(1nnnnnnyxfhyhxhfy
20、y31 的R-K公式(3.6)的局部误差不可能提高到 .2r)(4hO 把 多展开一项，从（3.8）的 看到展开式中 的项是不能通过选择参数消掉的.2Kny ffffyxy 实际上要使 的项为零，需增加3个方程，要确定4个参数 ，这是不可能的.3h21221,cc 故 的显式R-K方法的阶只能是 ，而不能得到三阶公式.2r2p329.3.3 9.3.3 三阶与四阶显式三阶与四阶显式R-KR-K方法方法 要得到三阶显式R-K方法，必须 .此时（3.4），（3.5）的公式表示为 3r).,(),(),(),(232131321212213322111hKhKyhxfKhKyhxfKyxfKKcKc
21、Kchyynnnnnnnn（3.11）其中 及 均为待定参数，公式（3.11）的局部截断误差为 321,ccc32313212,.)()(33221111KcKcKchxyxyTnnn只要将 按二元函数泰勒展开，使 ，可32,KK)(41hOTn33得待定参数满足方程.61,31,21,13223233222332232313212321cccccccc（3.12）这是8个未知数6个方程的方程组，解也不是唯一的.可以得到很多公式.34满足条件（3.12）的公式（3.11）统称为三阶R-K公式.一个常见的公式为).2,(),2,2(),(),4(62131213211hKhKyhxfKKhyhx
22、fKyxfKKKKhyynnnnnnnn此公式称为库塔库塔三阶方法.继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个：35).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn 四阶龙格-库塔方法的每一步需要计算四次函数值 ，可以证明其局部截断误差为 5()O h(,)f x y36 例例3 3 设取步长 ，从 直到 用四阶龙格-库塔方法求解初值问题 2.0h0 x1x.1)0(),10(2yxyxyy 解解 这里，经典的四阶龙格-库塔公式（3
23、.13）具有形式 37.)(2,222,222,2),22(6334223112143211hKyhxhKyKKhyhxKhyKKhyhxKhyKyxyKKKKKhyynnnnnnnnnnnnnn387321.17321.10.16125.16125.18.04832.14833.16.03416.13417.14.01832.11832.12.0)(39nnnxyyx计算结果表 比较例3和例2的计算结果，显然以龙格-库塔方法的精度为高.虽然四阶龙格-库塔方法的计算量（每一步要4次计算函数 ）比改进的欧拉方法（它是一种二阶龙格-库塔方法，每一步只要2次计算函数 ）大一倍，但由于这里放大了步长
24、，表9-3和表9-2 所耗费的计算量几乎相同.ff)2.0(h 龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质.反之，如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法.399.3.4 9.3.4 变步长的龙格变步长的龙格-库塔方法库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了.步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累.因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题.在选择步长时，需要考虑两个问题：1 怎样衡量和检验计算结果的
25、精度？2 如何依据所获得的精度处理步长？考察经典的四阶龙格-库塔公式（3.13）.40 从节点 出发，先以 为步长求出一个近似值，记为 ，由于公式的局部截断误差为 ，故有 nxh)(1hny)(5hO,)(5)(11chyxyhnn（3.14）然后将步长折半，即取 为步长从 跨两步到 ，再求得一个近似值 ，每跨一步的截断误差是 ，因此有 2hnx1nx)2(1hny52 hc,22)(5)2(11hcyxyhnn（3.15）比较（3.14）式和（3.15）式我们看到，步长折半后，误41差大约减少到 ，即有161.161)()()(11)2(11hnnhnnyxyyxy由此易得下列事后估计式.151)()(1)2(1)2(11hnhnhnnyyyxy这样，可以通过检查步长，折半前后两次计算结果的偏差)(1)2(1hnhnyy42来判定所选的步长是否合适，具体地说，将区分以下两种情况处理：1.对于给定的精度 ，如果 ，反复将步长折半进行计算，直至 为止，这时取最终得到的 作为结果；)2(1hny 2.如果 ，反复将步长加倍，直到 为止，这时再将步长折半一次，就得到所要的结果.这种通过加倍或折半处理步长的方法称为变步长方法变步长方法.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.
