第8章常微分方程初值问题数值解法8.1 引言8.2 欧拉方法8.3 龙格-库塔方法8.4 单步法的收敛性与稳定性8.5 线性多步法8.1 引 言考虑一阶常微分方程的初值问题00(,),[,],().y f x y x a b y x y '=∈=（1.1） （1.2）如果存在实数 ，使得121212(,)(,).,R f x y f x y L y y y y -≤-∀∈（1.3）则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点)(x y<<<<<+121n n x x x x 上的近似值 .,,,,,121+n n y y y y 相邻两个节点的间距 称为步长.n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ),2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”.即求解过程顺着节点排列的次序一步一步地向前推进.00(,),[,],().y f x y x a b y x y '=∈=描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式.1+n y 一类是计算时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值， 1+n y k 11,,,+--k n n n y y y 称为 步法.k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题.n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式.),(y x f y ='8.2 简单的数值方法8.2.1 欧拉法与后退欧拉法积分曲线上一点 的切线斜率等于函数 的 值.),(y x ),(y x f 如果按函数 在平面上建立一个方向场，那么， ),(y x f xy 积分曲线上每一点的切线方向均与方向场在该点的方向相一致.在 平面上，微分方程的解 称 作它的积分曲线.xy )(x y y =),(y x f y ='基于上述几何解释，从初始点 出发， ),(000y x P 先依切线 在该点的方向推进到 上一点 ，然后再 1x x =1P 从 依切线 的方向推进到 上一点 ，循此前进 1P 2x x =2P 做出一条折线 （图8-1).210P P P 图8-179一般地，设已做出该折线的顶点 ，过 依 切线 的方向再推进到 ，显然两个顶点的坐标有关系 ),(n n n y x P nP ),(111+++n n n y x P 1,+n n P P 11n n n n y y x x ++-=-反解，得).,(1n n n n y x hf y y +=+（2.1）这就是著名的欧拉（Euler ）方法. 欧拉法实际上是对常微分方程中的导数用差商近似，即1()()n n y x y x h+-≈直接得到的.(,)d d n n x y y x =(,)n n f x y 00(,),[,],().y f x y x a b y x y '=∈=1n n y y h +-=()(,())n n n y x f x y x '=111(,)n n n n y y hf x y +=+例欧拉方法13若初值已知，则依公式（2.1）可逐步算出 0y ),,(0001y x hf y y +=),,(1112y x hf y y +=).,(1n n n n y x hf y y +=+（2.1）例1 求解初值问题⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）解 欧拉公式的具体形式为).2(1nnn n n y x y h y y -+=+取步长 ，计算结果见表8-1.1.0=h 初值问题（2.2）的解为 ，按这个解析式 子算出的准确值 同近似值 一起列在表8-1中，两者相比较可以看出欧拉方法的精度不高.x y 21+=)(n x y n y 1(,)n n n n y y hf x y +=+12y x=+81()()0.1 1.1000 1.09540.6 1.5090 1.48320.2 1.1918 1.18320.7 1.5803 1.54920.3 1.2774 1.26490.8 1.6498 1.61250.4 1.3582 1.34160.9 1.7178 1.67330.51.43511.41421.01.78481.7321n nn n nn x y y x x y y x -表计算结果对比 还可以通过几何直观来考察欧拉方法的精度.假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图8-2）. )(n n x y y =nP )(x y y =1+n n P P )(x y y =nP图8-2从图形上看，这样定出的顶点 显著地偏离了原来 的积分曲线，可见欧拉方法是相当粗糙的.1+n P 为了分析计算公式的精度，通常可用泰勒展开将)(1+n x y在处展开，则有n x)()(1h x y x y n n +=+在 的前提下， )(n n x y y =)())(,(),(n n n n n x y x y x f y x f '==11()n n y x y ++-=称为此方法的局部截断误差.()()n n y x y x h '=++于是可得欧拉法（2.1）的误差 （2.3）),(22n x y h ''≈如果对方程 从 到积分，得 n x 1+n x ),(y x f y =').,(1n n n n y x hf y y +=+（2.1） 2()2n h y ξ''1(,)n n n x x ξ+∈1(,)n n n n y y hf x y +=+(上一值精确相等时，估计下一值误差!) 2()2n h y ξ''1()n y x +=（2.4）1()(,())n nx n xy x f t y t dt++⎰18右端积分用左矩形公式 近似.))(,(n n x y x f h 再以 代替 n y ),(n x y 如果在（2.4）中右端积分用右矩形公式 ))(,(11++n n x y x f h 111(,)n n n n y y hf x y +++=+（2.5）称为后退的欧拉法.代替也得到（2.1）， )(1+n x y 1+n y 局部截断误差也是（2.3）. 近似，则得另一个公式它也可以通过利用均差近似导数 ，即11()()n n n ny x y x x x ++-≈-)(1+'n x y 1(,)n n n n y y hf x y +=+2()2n hy ξ''1()(,())n nx n x y x f t y t dt++⎰111()(,())n n n y x f x y x +++'=欧拉公式是关于 的一个直接的计算公式，这类公式称作是显式的；1+n y 后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的.1+n y 1+n y 1(,)n n n n y y hf x y +=+111(,)n n n n y y hf x y +++=+隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化. 设用欧拉公式),()0(1n n n n y x f h y y+=+给出迭代初值 ，用它代入（2.5）式的右端，使之转化 为显式，直接计算得)0(1+n y ),,()0(11)1(1++++=n n n n yx f h y y然后再用 代入（2.5）式，又有)1(1+n y ).,()1(11)2(1++++=n n n n yx f h y y111(,)n n n n y y hf x y +++=+（2.5）21如此反复进行，得).,1,0(),,()(11)1(1=+=++++k yx f h y yk n n n k n （2.6）由于 对 满足李普希茨条件(1.3). ),(y x f y 由(2.6)减(2.5)得(1)11k n n yy +++-=.1)(1++-≤n k n y y hL 由此可知，只要迭代法(2.6)就收敛到解 . 1<hL 1+n y 从积分公式可以看到后退欧拉方法的公式误差与欧拉法是相似的.111(,)n n n n y y hf x y +++=+（2.5） ()1111(,)(,)k n n n n h f x y f x y ++++-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈8.2.2 梯形方法 若用梯形求积公式近似等式（2.4）右端的积分并分别用 代替 则可得到比欧拉法 精度高的计算公式1,+n n y y ),(),(1+n n x y x y 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） 称为梯形方法. 梯形方法是隐式单步法，可用迭代法求解.⎰+1))(,(n nx x dtt y t f .))(,()()(11⎰++=+n nx xn n dt t y t f x y x y （2.4） 1()(,())n nx n xy x f t y t dt++⎰111(,())[(,)(,)]2n n x n n n n x hf t y t dt f x y f x y +++≈+⎰23⎪⎪⎩⎪⎪⎨⎧+=+);,()0(1n n n ny x f h y y 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得)],(),([2)(11)1(1k nn n n n k n y x f y x f h y y ++++++=（2.8） ).,2,1,0( =k 同后退的欧拉方法一样，仍用欧拉方法提供迭代初值， 则梯形法的迭代公式为111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-如果选取 充分小，使得h 12hL<则当 时有 ， ∞→k 1)(1++→n k n y y 此时迭代过程是收敛的.于是有(1)11k n n y y+++-≤式中 为关于 的李普希茨常数. ),(y x f y L (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈()112k n n hL y y ++-计算方法及MATLAB 实现8.2.3 改进欧拉公式梯形方法虽然提高了精度，但其算法复杂. 在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.),(y x f 为了控制计算量，通常只迭代一两次就转入下一步的 计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 , 1n y + 而迭代又要反复进行若干次，计算量很大，而且往往难以预测.称之为预测值，)],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=计算方法及MATLAB 实现这样建立的预测-校正系统通常称为改进的欧拉公式： 预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正值. 1+n y 1+n y 预测 ),,(1n n n n y x f h y y +=+校正)].,(),([2111+++++=n n n n n n y x f y x f hy y （2.9）也可以表为下列平均化形式),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(1y y y +=111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） )],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=(参阅教材392页!)计算方法及MATLAB 实现例2 用改进的欧拉方法求解初值问题（2.2）.解 这里 改进的欧拉公式为⎪⎪⎪⎩⎪⎪⎪⎨⎧-+=),2(nnn n p y x y h y y ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）2(,)(),xf x y y y=-),2(1p n p n c y x y h y y +-+=).(211c p n y y y +=+),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+计算方法及MATLAB 实现282(,)xf x y y y=-),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+取 ，计算结果见表8-2.1.0=h 同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.()()0.11.0959 1.09540.6 1.4860 1.48320.2 1.1841 1.18320.7 1.5525 1.54920.3 1.2662 1.26490.8 1.6153 1.61250.4 1.3434 1.34160.9 1.6782 1.67330.51.41641.41421.01.73791.7321n n n n n n x y y x x y y x -表8212y x=+9.2.4 单步法的局部截断误差与阶 初值问题(1.1)，(1.2)的单步法可用一般形式表示为),,,,(11h y y x h y y n n n n n +++=ϕ（2.10）其中多元函数 与 有关，ϕ),(y x f 当 含有 时，方法是隐式的，若不含 则为显式方法，ϕ1+n y 1+n y ),,,(1h y x h y y n n n n ϕ+=+（2.11） 称为增量函数， ),,(h y x ϕ).,(),,(y x f h y x =ϕ所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ).,(1n n n n y x hf y y +=+（2.1） )(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈对一般显式单步法则可如下定义. 定义1 设 是初值问题(1.1)，(1.2)的准确解， 称 )(x y y =)),(,()()(11h x y x h x y x y T n n n n n ϕ--=++（2.12） 为显式单步法（2.11）的局部截断误差. 之所以称为局部的，是假设在 前各步没有误差.1+n T n x 当 时，计算一步，则有 )(n n x y y =11()n n y x y ++-=1()()(,(),)n n n n y x y x h x y x h ϕ+=--⎩⎨⎧=='.)(),,(00y x y y x f y （1.1）（1.2） ),,,(1h y x h y y n n n n ϕ+=+（2.11）1()[(,,)]n n n n y x y h x y h ϕ+-+1.n T +=1(,)n n n n y y hf x y +=+).,(),,(y x f h y x =ϕ在前一步精确的情况下用显式单步公式计算产生的公式误差. 根据定义，欧拉法的局部截断误差 ))(,()()(11n n n n n x y x hf x y x y T --=++即为（2.3）的结果. 这里 称为局部截断误差主项. )(22n x y h '')(2h O T =局部截断误差可理解为用显式单步方法计算一步的误差，即 ()()n n y x h y x =+-22()()2n h y x o h ''=+显然 ),,,(1h y x h y y n n n n ϕ+=+（2.11）)(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈()n hy x '-(,)y f x y '=1()()n n y x y x h +=+(泰勒公式) (小o 高阶，大O 同阶)计算方法及MATLAB 实现 定义2 设 是初值问题(1.1)，(1.2)的准确解， 若存在最大整数使显式单步法(2.11)的局部截断误差满足 )(x y p ),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ（2.13） 则称显式单步法具有 阶精度. p 若将（2.13）展开式写成),())(,(211++++=p p n n n h O h x y x T ψ则 称为局部截断误差主项. 1))(,(+p n n hx y x ψ 以上定义对隐式单步法（2.10）也是适用的.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ),,,,(11h y y x h y y n n n n n +++=ϕ（2.10） 1(,,)n n n n y y h x y h ϕ+=+(小o 高阶，大O 同阶)对后退欧拉法（2.5）其局部截断误差为))(,()()(1111++++--=n n n n n x y x hf x y x y T 这里 ，是1阶方法，局部截断误差主项为 . 1=p )(22n x y h ''-23()()()2n n h hy x y x O h '''=++).()(232h O x y h n +''-=),,(111++++=n n n n y x hf y y （2.5）2[()()()]n n h y x hy x O h '''-++对梯形法（2.7）有)]()([2)()(111+++'+'--=n n n n n x y x y h x y x y T 所以梯形方法是二阶的，其局部误差主项为 )(123n x y h '''-23()()()23!n n n h h hy x y x y x ''''''=++).()(1243h O x y h n +'''-=)],,(),([2111+++++=n n n n n n y x f y x f h y y （2.7） 24[()()()()]()22n n n n h h y x y x hy x y x O h '''''''-++++3364h h -8.3龙格-库塔方法8.3.1 显式龙格-库塔法的一般形式 上节给出了显式单步法的表达式),,,(1h y x h y y n n n n ϕ+=+其局部截断误差为),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ对欧拉法 ，即方法为阶， )(21h O T n =+1=p ))].,(,(),([21n n n n n n n n y x hf y h x f y x f h y y ++++=+（3.1）若用改进欧拉法，它可表为此时增量函数 ))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ（3.2）与欧拉法的 相比，增加了计算一个 右函数 的值，可望 .),(),,(n n n n y x f h y x =ϕf 2=p 若要使得到的公式阶数 更大， 就必须包含更多的值.p ϕf ,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）从方程 等价的积分形式（2.4），即),(y x f y ='若要使公式阶数提高，就必须使右端积分的数值求积公式 精度提高，必然要增加求积节点.为此可将（3.3）的右端用求积公式表示为线性组合1(,())n n x x f x y x dx +≈⎰点数 越多，精度越高， r 上式右端相当于增量函数 ， ),,(h y x ϕ为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为),,,(1h y x h y y n n n n ϕ+=+（3.4）其中1(,())r i n i n i i h c f x h y x h λλ=++∑,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）(可变步长因子!) (组合系数!),),,(1∑==r i i i n n K c h y x ϕ),,(1n n y x f K =),(11∑-=++=i j j ij n i n i K h y h x f K μλ（3.5）,,,2r i =这里 均为常数. ij i i c μλ,, (3.4)和(3.5)称为 级显式龙格-库塔（Runge-Kutta ）法， r 简称R-K 方法.当 时，就是欧拉法，),(),,(,1n n n n y x f h y x r ==ϕ此时方法的阶为. 1=p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，2=r ),,,(1h y x h y y n n n n ϕ+=+（3.4） ),,(1h y x h y y n n n n ϕ+=+))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ).,(1n n n n y x hf y y+=+1(,,)n n n n y y h x y h ϕ+=+(步长调节因子) (组合系数) (斜率组合系数)8.3.2 二阶显式R-K 方法 1(,,)rn n i ii x y h c K ϕ==∑1(,,)n n n n y y h x y h ϕ+=+ 对 的R-K 方法，计算公式如下2=r 11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩（3.6）这里 均为待定常数.21221,,,μλc c ),(11∑-=++=i j j ij n i n i K h y h x f K μλ若取 得计算公式,2/1,1,021221====μλc c 12n n y y hK +=+称为中点公式，相当于数值积分的中矩形公式.上式也可表示为1(,(,))22n n n n n n h hy y hf x y f x y +=+++1(,)n n K f x y =（3.10）21(,)22n n h hK f x y K =++)).,(2,2(1n n n n n n y x f hy h x hf y y +++=+11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩继续上述过程，经过较复杂的数学演算，可以导出各 种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为 . )(5h O 四阶龙格-库塔方法的每一步需要计算四次函数值 ，f ⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+).,(),2,2(),2,2(),,(),22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n n n n n n n （3.13） (凸线性组合!) (中点公式!) (欧拉公式!)例 设取步长 ，从 到 用四阶龙格-库塔方法求解 初值问题 ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y xy y 解 这里 ，经典的四阶龙格-库塔公式为 yx y y x f 2),(-=),22(643211K K K K hy y n n ++++=+,2nx y K -=2.0=h 0=x 1=x 112341213243(22),6(,),(,),22(,),22(,).n n n nn n n n n n h y y K K K K K f x y h h K f x y K h hK f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩计算方法及MATLAB 实现,222223K h y h x K h y K n n n ++-+=,222112K h y h x K hy K n n n ++-+=.)(2334hK y h x hK y K n n n ++-+= 表8-3列出了计算结果，同时了出了相应的精确解. 龙格-库塔方法的精度较高.112341213243(22)6(,),(,),22(,),22(,).n n n nn n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩yxy y x f 2),(-=,21nnn y x y K -=83()0.2 1.1832 1.18320.4 1.3417 1.34160.6 1.4833 1.48320.8 1.6125 1.61251.01.73211.7321n n n x y y x 表计算结果。