第9章 常微分方程初值问题数值解法

1、上页上页上页上页上页上页下页下页下页下页下页下页第第9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法 9.1 引言引言 9.2 简单的数值方法简单的数值方法 9.3 龙格龙格- -库塔方法库塔方法 9.4 单步法的收敛性与稳定性单步法的收敛性与稳定性 9.5 线性多步法线性多步法 9.6 线性多步法的收敛性与稳定性线性多步法的收敛性与稳定性 9.7 一阶方程组与刚性方程组一阶方程组与刚性方程组上页上页上页上页上页上页下页下页下页下页下页下页9.1 引引 言言 科学技术中很多问题都可用微分方程的科学技术中很多问题都可用微分方程的定解问定解问题题来描述，主要有来描述，主要有初值问题初值

2、问题与与边值问题边值问题两大类，两大类，本章本章只考虑初值问题只考虑初值问题. 常微分方程初值问题中最简单的例常微分方程初值问题中最简单的例子是人口模型，设某特定区域在子是人口模型，设某特定区域在t0时刻人口为时刻人口为y(t0)=y0已知的，该区域的人口自然增长率为已知的，该区域的人口自然增长率为 ，人口增长与，人口增长与人口总数成正比，所以人口总数成正比，所以t时刻的人口总数时刻的人口总数y(t)满足以下满足以下微分方程微分方程 .)(),(00ytytyy 上页上页上页上页上页上页下页下页下页下页下页下页很多物理系统与时间有关，从卫星运行轨道到单很多物理系统与时间有关，从卫星运行轨道到单

3、摆运动，从化学反应到物种竞争都是随时间的延续而摆运动，从化学反应到物种竞争都是随时间的延续而不断变化的不断变化的. 解常微分方程是描述连续变化的数学语解常微分方程是描述连续变化的数学语言，微分方程的求解就是确定满足给定方程的可微函言，微分方程的求解就是确定满足给定方程的可微函数数y(t)，研究它的数值方法是本章的主要目的，研究它的数值方法是本章的主要目的. 考虑考虑一一阶常微分方程的初值问题阶常微分方程的初值问题 )2 . 1(.)()1 . 1(,),(000yxybxxyxfy则称则称f关于关于y满足满足利普希茨利普希茨(Lipschitz)条件条件，L称为称为y的的利普希茨常数利普希茨常

4、数(简称简称Lips.常数常数).如果存在实数如果存在实数L0，使得，使得)3 . 1(.,),(),(212121RyyyyLyxfyxf 上页上页上页上页上页上页下页下页下页下页下页下页 定理定理1 设设f在区域在区域D=(x, y)|a x b,y R上连续上连续, 关关于于y满足利普希茨条件，则对任意满足利普希茨条件，则对任意x0 a, b, y0 R，常，常微分方程初值问题微分方程初值问题(1.1)式和式和(1.2)式当式当x a, b时存在时存在唯一的连续可微解唯一的连续可微解y(x) .解的存在唯一性定理是常微分方程理论的基本内解的存在唯一性定理是常微分方程理论的基本内容，也是数

5、值方法的出发点，此外还要考虑方程的解容，也是数值方法的出发点，此外还要考虑方程的解对扰动的敏感性，它有以下结论对扰动的敏感性，它有以下结论. 定理定理2 设设f在区域在区域D (如定理如定理1所定义所定义) 上连续上连续, 且且关于关于y满足利普希茨条件，设初值问题满足利普希茨条件，设初值问题.)(),()(0sxyyxfxy 的解为的解为y(x, s)，则，则.),(),(21210ssesxysxyxxL 上页上页上页上页上页上页下页下页下页下页下页下页这个定理表明解对初值依赖的敏感性，它与右端这个定理表明解对初值依赖的敏感性，它与右端函数函数f有关，当有关，当f的的Lips.常数常数L比

6、较小时，解对初值和比较小时，解对初值和右端函数相对不敏感，可视为好条件右端函数相对不敏感，可视为好条件. 若若L较大则可较大则可认为坏条件，即病态问题认为坏条件，即病态问题.如果右端函数可导，由中值定理有如果右端函数可导，由中值定理有.,),(),(),(212121之间之间在在yyyyyxfyxfyxf 若假定若假定 在域在域D内有界内有界, 设设 , 则则yyxf ),(Lyyxf ),(.),(),(2121yyLyxfyxf 它表明它表明f满足利普希茨条件，且满足利普希茨条件，且L的大小反映了右端函的大小反映了右端函数数f关于关于y变化的快慢，刻画了初值问变化的快慢，刻画了初值问(1.

7、1)式和式和(1.2)式式是否为好条件是否为好条件. 这在数值求解中也是很重要的这在数值求解中也是很重要的.上页上页上页上页上页上页下页下页下页下页下页下页 虽然求解常微分方程有各种各样的解析方法，但虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法题中归结出来的微分方程主要靠数值解法. 所谓所谓数值解法数值解法, 就是寻求解就是寻求解y(x)在一系列离散节点在一系列离散节点 121nnxxxx上的近似值上的近似值 y1,y2,yn,yn+1,. 相邻两个节点的间距相邻两个节

8、点的间距hn=xn+1- -xn称为称为步长步长. 今后如不特别说明，总是假定今后如不特别说明，总是假定 hi=h(i=1,2,)为为常数常数, 这时节点为这时节点为xn=x0+nh(i=0,1,2,) (等距节点等距节点).上页上页上页上页上页上页下页下页下页下页下页下页 初值问题的初值问题的数值解法数值解法有个有个基本特点基本特点，他们都采，他们都采取取“步进式步进式”，即求解过程顺着节点排列的次序一步，即求解过程顺着节点排列的次序一步一步地向前推进一步地向前推进. 描述这类算法，只要给出用已知信描述这类算法，只要给出用已知信息息yn,yn- -1,yn- -2,计算计算yn+1的递推公式

9、的递推公式. 本章首先要对常微分方程本章首先要对常微分方程(1.1)离散化，建立求解离散化，建立求解数值解的递推公式数值解的递推公式. 一类是计算一类是计算yn+1时只用到前一点的时只用到前一点的值值yn，称为，称为单步法单步法. 另一类是用到另一类是用到yn+1前面前面 k 点的值点的值yn,yn-1, yn-k+1，称为，称为k步法步法. 其次，要研究公式的其次，要研究公式的局局部截断误差部截断误差和和阶阶，数值解，数值解yn与精确解与精确解y(xn)的的误差估计误差估计及及收敛性收敛性，还有递推公式的，还有递推公式的计算稳定性计算稳定性等问题等问题.上页上页上页上页上页上页下页下页下页下

10、页下页下页9.2 简单的数值方法简单的数值方法9.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法 我们知道，在我们知道，在xy平面上，微分方程平面上，微分方程(1.1)式式的解的解y=f(x)称作它的称作它的积分曲线积分曲线，积分曲线积分曲线上一点上一点(x, y)的切的切线斜率等于函数线斜率等于函数f(x, y)的值的值. 如果按如果按f(x, y)在在xy平面上平面上建立一个方向场，那么，建立一个方向场，那么，积分曲线积分曲线上每一点的切线上每一点的切线方向均与方向场在该点的方向相一致方向均与方向场在该点的方向相一致.基于上述几何解释，我们从初始点基于上述几何解释，我们从初始点P0(x0,

11、y0)出发出发,先依方向场在该点的方向推进到先依方向场在该点的方向推进到x=x1上一点上一点P1，然后，然后再从再从P1点依方向场在该点的方向推进到点依方向场在该点的方向推进到 x=x2 上一点上一点P2 , 循环前进做出一条循环前进做出一条折线折线P0 P1 P2.上页上页上页上页上页上页下页下页下页下页下页下页 一般地，设已做出该折线的顶点一般地，设已做出该折线的顶点Pn, ,过过Pn(xn, yn)依方向场的方向再推进到依方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个，显然两个顶点顶点Pn, ,Pn+1的坐标有关系的坐标有关系),(,()(),(111nnnnnnnnnn

12、nxyxfxyyxfhyyxxyy ) 1 .2 (),(1nnnnyxhfyy 这就是著名的这就是著名的( (显式显式) )欧拉欧拉( (Euler) )公式公式. . 若初值若初值y0已已知，则依公式知，则依公式(2.1)可逐次逐步算出各点数值解可逐次逐步算出各点数值解. .即即),(0001yxhfyy ),(1112yxhfyy 上页上页上页上页上页上页下页下页下页下页下页下页 例例1 用欧拉公式求解初值问题用欧拉公式求解初值问题)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步长取步长h=0.1，欧拉公式的具体形式为，欧拉公式的具体形式为)2(1nnnnnyxyhyy

13、其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy上页上页上页上页上页上页下页下页下页下页下页下页依次计算下去，依次计算下去，部分计算结果部分计算结果见下表见下表. xy21 与准确解与准确解 相比，可看出欧拉公式的计算结相比，可看出欧拉公式的计算结果精度很差果精度很差. xn 欧拉公式数值解欧拉公式数值解yn准确解准确解y(xn) 误差误差 0.2 0.4 0.6 0.8 1.0 1.1918

14、18 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719上页上页上页上页上页上页下页下页下页下页下页下页 欧拉公式具有明显的几何意义欧拉公式具有明显的几何意义, , 就是就是用折线近似用折线近似代替方程的解曲线代替方程的解曲线，因而常称公式，因而常称公式(2.1)为为欧拉折线法欧拉折线法. .( )yy xxynx1nxnp1np1np x 还可以通过几何直观来考察欧拉方法的精度还可以通过几何直观

15、来考察欧拉方法的精度. .假假设设yn=y(xn), ,即顶点即顶点Pn落在积分曲线落在积分曲线y=y(x)上，那么，上，那么，按欧拉方法做出的折线按欧拉方法做出的折线PnPn+1便是便是y=y(x)过点过点Pn的切线的切线. .从图形上看从图形上看, ,这这样定出的顶点样定出的顶点Pn+1显著显著地偏离了原来的积分曲地偏离了原来的积分曲线，可见欧拉方法是线，可见欧拉方法是相相当粗糙当粗糙的的. .上页上页上页上页上页上页下页下页下页下页下页下页 为了分析计算公式的精度，通常可用泰勒展开为了分析计算公式的精度，通常可用泰勒展开将将y(xn+1)在在xn处展开，则有处展开，则有).,()(2),

16、()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y(xn)的前提下，的前提下，f(xn,yn )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得欧拉法可得欧拉法(2.1)的的公式误差公式误差为为) 3 . 2(),(2)(2)(2211nnnnxyhyhyxy 称为此方法的称为此方法的局部截断误差局部截断误差. .上页上页上页上页上页上页下页下页下页下页下页下页 如果对方程如果对方程(1.1)从从xn到到xn+1积分，得积分，得) 4 . 2(.)(,()()(11 nnxxnndttytfxyxy右端积分

17、用右端积分用左矩形公式左矩形公式hf(xn,y(xn)近似，再以近似，再以yn代替代替y(xn)，yn+1代替代替y(xn+1)也得到欧拉公式也得到欧拉公式(2.1)，局部截，局部截断误差也是断误差也是(2.3).称为称为( (隐式隐式) )后退的欧拉公式后退的欧拉公式. .它也可以通过利用均差它也可以通过利用均差近似导数近似导数y (xn+1)，即，即 如果右端积分用如果右端积分用右矩形公式右矩形公式hf(xn+1,y(xn+1)近似，近似，则得到另一个公式则得到另一个公式) 5 . 2 (),(111 nnnnyxhfyy直接得到直接得到. .),(,()()()(1111 nnnnnnn

18、xyxfxyxxxyxy上页上页上页上页上页上页下页下页下页下页下页下页 后退的欧拉公式与欧拉公式有着本质的区别后退的欧拉公式与欧拉公式有着本质的区别, 后后者是关于者是关于yn+1的一个直接计算公式，这类公式称作是的一个直接计算公式，这类公式称作是显式的显式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它实际上，它实际上是关于是关于yn+1的一个函数方程的一个函数方程, ,这类方程称作是这类方程称作是隐式的隐式的. . 显式显式与与隐式隐式两类方法各有特点，考虑到数值稳两类方法各有特点，考虑到数值稳定性等其他因素，人们有时需要选用定性等其他因素，人们有时需要选用隐式隐式方法

19、，但方法，但使用使用显式显式算法远比算法远比隐式隐式方便方便. . 隐式方程隐式方程(2.5)通常用迭代法求解，而迭代过程通常用迭代法求解，而迭代过程的实质是的实质是逐步逐步显式化显式化. .上页上页上页上页上页上页下页下页下页下页下页下页 设用欧拉公式设用欧拉公式),() 0(1nnnnyxhfyy 给出迭代初值给出迭代初值 ，用它代入，用它代入(2.5)式的式的右端，使之转右端，使之转化为显式，直接计算得化为显式，直接计算得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(

20、1 nnnnyxhfyy如此反复进行，得如此反复进行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn上页上页上页上页上页上页下页下页下页下页下页下页由于由于f(x, y)对对y满足满足Lipschitz条件条件(1.3). 由由(2.6)减减(2.5)得得),(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL , ,我们反复将步长我们反复将步长折半计算折半计算, ,直至直至 为为止止, ,这时取最终得到的这时取最终得到的 作为作为结果；结果； 21hny(2) 如果如果 为

21、止，这时再将步长折半计算一次，就得到所为止，这时再将步长折半计算一次，就得到所要的结果要的结果. .变步长方法还可利用变步长方法还可利用p阶与阶与p+1阶公式的局部截断阶公式的局部截断误差得到误差控制与变步长的具体方法误差得到误差控制与变步长的具体方法, ,参见文献参见文献8.8.上页上页上页上页上页上页下页下页下页下页下页下页9.4 单步法的收敛性与稳定性单步法的收敛性与稳定性9.4.1 收敛性与相容性收敛性与相容性 数值解法的基本思想是，通过某种离散化手段数值解法的基本思想是，通过某种离散化手段将微分方程将微分方程(1.1)转化为差分方程转化为差分方程,如单步法如单步法(2.11),即即

22、) 1 . 4(. ),(1hyxhyynnnn 它在点它在点xn处的解为处的解为yn，而初值问题，而初值问题(1.1), (1.2) 在点在点xn处的精确解为处的精确解为y( (xn) )，记，记en=y(xn)- -yn称为称为整体截断误整体截断误差差. 收敛性收敛性就是讨论当就是讨论当 x=xn 固定且固定且 时时en0的问题的问题.00 nxxhn上页上页上页上页上页上页下页下页下页下页下页下页 定义定义3 若一种数值方法对于固定的若一种数值方法对于固定的xn=x0+nh, 当当h0时有时有yny(xn)，其中，其中y(x)是是(1.1),(1.2)的准确解，的准确解，则称该方法是则称

23、该方法是收敛收敛的的. 显然数值方法收敛是指显然数值方法收敛是指en=y(xn)- -yn0，对单步，对单步法法(4.1)有下述收敛性定理：有下述收敛性定理： 定理定理3 假设单步法假设单步法(4.1)具有具有p阶精度，且增量函阶精度，且增量函数数 (x, y ,h)关于关于y满足利普希次条件满足利普希次条件又设初值又设初值y0是准确的是准确的, 即即y0=f(x0), 则其则其整体截断误差整体截断误差)2 . 4(.),(),(yyLhyxhyx ) 3 . 4().()(pnnhOyxy 上页上页上页上页上页上页下页下页下页下页下页下页 证明证明 设以设以yn+1表示取表示取yn=y(xn

24、)用公式用公式(4.1)求得的求得的结果，即结果，即) 4 . 4(),),(,()(1hxyxhxyynnnn 则则y(xn)- -yn+1为局部截断误差，由于所给方法具为局部截断误差，由于所给方法具有有p阶精度，按定义阶精度，按定义2，存在定数，存在定数C ，使，使.)(111 pnnChyxy又由式又由式(4.1)与与(4.4)，得，得. ),(),(,()(11hyxhxyxhyxyyynnnnnnnn 上页上页上页上页上页上页下页下页下页下页下页下页利用利普希次条件利用利普希次条件(4.2)，有，有.)()1 (11nnnnyxyhLyy 从而有从而有.)()1 ()()(11111

25、11 pnnnnnnnnChyxyhLyxyyyyxy 即对整体截断误差即对整体截断误差en=y(xn)- -yn成立下列成立下列递推关系式递推关系式) 5 . 4(.)1 (11 pnnChehLe 据此不等式反复递推，可得据此不等式反复递推，可得) 6 . 4(.1)1()1 (0 npnnhLLChehLe 上页上页上页上页上页上页下页下页下页下页下页下页由此可以断定，如果初值是准确的，即由此可以断定，如果初值是准确的，即e0=0，则，则(4.3)式成立式成立. 定理证毕定理证毕.再注意到当再注意到当x=x0+nh T时时 ,)1 ( TLnhLneehL 最终得下列估计式最终得下列估计

26、式) 7 . 4().1(0 TLpTLneLCheee依据这一定理，判断单步法依据这一定理，判断单步法(4.1)的收敛性，归结的收敛性，归结为验证增量函数为验证增量函数 能否满足利普希次条件能否满足利普希次条件(4.2).对于欧拉方法，由于其增量函数对于欧拉方法，由于其增量函数 就是就是f(x, y), 故当故当f(x, y)关于关于y满足利普希次条件时它是收敛的满足利普希次条件时它是收敛的.上页上页上页上页上页上页下页下页下页下页下页下页再考察改进的欧拉方法，其增量函数再考察改进的欧拉方法，其增量函数 已由已由(3.2)式给出，假定式给出，假定f(x, y)关于关于y满足利普希次条件，即满

27、足利普希次条件，即 ),(),(hyxhyx 这时有这时有.),(),(yyLyxfyxf ),(),(21yxfyxf .2121),(),(212121),(),(2121),(,(),(,(212yyLhLyyhLyyLyxfyxfhLyyLyyLyxhfyyxhfyLyyLyxhfyhxfyxhfyhxf 上页上页上页上页上页上页下页下页下页下页下页下页即即.21),(),(yyLhLhyxhyx 设限定步长设限定步长hh0(h0为定数为定数)，上式表明，上式表明关于关于y的利普的利普希次常数为希次常数为.21 LhLL 因此改进的欧拉方法也是收敛的因此改进的欧拉方法也是收敛的.类似地

28、类似地, 不难验证其它龙格不难验证其它龙格- -库塔方法的收敛性库塔方法的收敛性.上页上页上页上页上页上页下页下页下页下页下页下页 定理定理3表明表明p 1时单步法收敛时单步法收敛, 并且当并且当y(x)是初值是初值问题问题(1.1),(1.2)的解的解, (4.1)具有具有p阶精度时阶精度时, 则有展开式则有展开式),(,()()(1hxyxhxyhxyTn ) 0),(,() 0),(,(2)()(2 hxyxxyxhhxyhxyx ).()0),(,()(2hOxyxxyh 所以所以p 1的充分必要条件是的充分必要条件是 ，而而 ，于是可给出如下定义：，于是可给出如下定义：0) 0),(

29、,()( xyxxy )(,()(xyxfxy 上页上页上页上页上页上页下页下页下页下页下页下页 定义定义4 若单步法若单步法(4.1)的增量函数的增量函数 满足满足 以上讨论表明以上讨论表明p阶方法阶方法(4.1)当当p 1时与时与(1.1), (1.2)相容，反之相容，反之相容方法至少是相容方法至少是1阶的阶的.),()0 ,(yxfyx 则称单步法则称单步法(4.1)与初值问题与初值问题(1.1),(1.2)相容相容. 相容性是指方法逼近微分方程相容性是指方法逼近微分方程(1.1),即即(1.1)离散离散化得到的数值方法化得到的数值方法, 当当h0时可得到时可得到y (x)=f(x, y

30、) .于是有下面定理于是有下面定理. 定理定理4 p阶方法阶方法(4.1)与初值问题与初值问题(1.1),(1.2)相容的相容的充分必要条件是充分必要条件是p 1. 于是由定理于是由定理3可知方法可知方法(4.1)收敛的收敛的充分必要条件充分必要条件是是此方法是相容此方法是相容的的.上页上页上页上页上页上页下页下页下页下页下页下页9.4.2 绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域 前面关于收敛性的讨论有个前提，必须假定数前面关于收敛性的讨论有个前提，必须假定数值方法本身的计算是准确的值方法本身的计算是准确的. 实际情形并不是这样，实际情形并不是这样，差分方程的求解还会有差分方程的求解还会

31、有计算误差计算误差. 譬如由于数字舍入譬如由于数字舍入而引起的而引起的小扰动小扰动. 这类小扰动在传播过程中会不会恶这类小扰动在传播过程中会不会恶性增长，以至于性增长，以至于“淹没淹没”了差分方程的了差分方程的“真解真解”呢？呢？这就是这就是差分方程的稳定性问题差分方程的稳定性问题. 在实际计算时，我们在实际计算时，我们希望某一步产生的扰动值，在后面的计算中希望某一步产生的扰动值，在后面的计算中能够被能够被控制控制，甚至是，甚至是逐步衰减逐步衰减的的.上页上页上页上页上页上页下页下页下页下页下页下页 定义定义5 若一种数值方法在节点值若一种数值方法在节点值yn上大小为上大小为 扰扰动，于以后各

32、节点值动，于以后各节点值ym(mn)上产生的偏差均不超过上产生的偏差均不超过 ，则称该方法是，则称该方法是稳定稳定的的. 下面以欧拉法为例考察计算稳定性下面以欧拉法为例考察计算稳定性. 例例4 用欧拉公式求解初值问题用欧拉公式求解初值问题 . 1)0(,100yyy 解解 用欧拉法解方程用欧拉法解方程y=- -100y 得得其准确解其准确解 是一个按指数曲线衰减很快的是一个按指数曲线衰减很快的函数函数.(见书见书p294的图的图9-3)xexy100)( .)1001(1nnyhy 上页上页上页上页上页上页下页下页下页下页下页下页若取步长若取步长h=0.025，则欧拉公式的具体形式为，则欧拉公

33、式的具体形式为.5 . 11nnyy 节点节点xn欧拉方法欧拉方法yn后退欧拉方法后退欧拉方法yn0.0250.0500.0750.100- -1.5 2.25 - -3.375 5.06250.28570.08160.02330.0067计算结果见表计算结果见表, , 明显计算过程不稳定明显计算过程不稳定, , 但取但取h=0.005, , yn+1=0 0.5yn, , 则计算过程稳定则计算过程稳定. .对后退的欧拉公式，取对后退的欧拉公式，取h=0.025时，则计算公式时，则计算公式为为yn+1=(1/3.5)yn . .计算结果见表计算结果见表, , 这时计算过程是稳这时计算过程是稳定

34、的定的. .上页上页上页上页上页上页下页下页下页下页下页下页 例题表明稳定性不但例题表明稳定性不但与方法有关与方法有关，也，也与步长与步长h有有关，当然关，当然与方程中的与方程中的f(x, y)有关有关. 为了只考察数值方为了只考察数值方法本身，通常只检验数值方法用于解法本身，通常只检验数值方法用于解模型方程模型方程的稳的稳定性，定性，模型方程模型方程为为)8 . 4(, yy 其中其中 为复数，这个方程分析较简单，对一般方程可为复数，这个方程分析较简单，对一般方程可以通过局部线性化化为这种形式，例如在以通过局部线性化化为这种形式，例如在(x, y)的邻的邻域，可展开为域，可展开为,)(,()

35、(,(),(),( yyyxfxxyxfyxfyxfyyx上页上页上页上页上页上页下页下页下页下页下页下页略去高阶项，再做变换即可得到略去高阶项，再做变换即可得到u = u的形式的形式. 对于对于m个方程的常微分方程组个方程的常微分方程组, 可线性化为可线性化为y =Ay, 这里矩这里矩阵阵A为为mm雅可比矩阵雅可比矩阵(fi/yj)，若矩阵，若矩阵A有有m个特征个特征值值 1, 2, m，其中，其中 i可能是复数，所以，为了使模可能是复数，所以，为了使模型方程结果能推广到常微分方程组，方程型方程结果能推广到常微分方程组，方程(4.8)式中式中 为复数为复数. 为保证微分方程本身的稳定性，还应

36、假定为保证微分方程本身的稳定性，还应假定Re( )0.下面先研究欧拉方法的稳定性下面先研究欧拉方法的稳定性. 模型方程模型方程y = y的欧拉公式为的欧拉公式为)9 . 4(.)1(1nnyhy 上页上页上页上页上页上页下页下页下页下页下页下页设在节点设在节点 yn 上有一扰动值上有一扰动值 n，它的传播使节点值，它的传播使节点值yn+1产生大小为的扰动值产生大小为的扰动值 n+1，假设用，假设用y*n=yn+ n，按欧拉，按欧拉公式得出公式得出 y*n+1=yn+1+ n+1 的计算过程不再有新的误差，的计算过程不再有新的误差，则扰动值满足则扰动值满足.)1(1nnh 可见扰动值满足原来的差

37、分方程可见扰动值满足原来的差分方程(4.9). 这样，如果差这样，如果差分方程的解是不增长的，即有分方程的解是不增长的，即有.1nnyy 则它就是稳定的则它就是稳定的. 这一论断对于下面将要研究的其它这一论断对于下面将要研究的其它方法同样适用方法同样适用.上页上页上页上页上页上页下页下页下页下页下页下页显然，为要保证差分方程显然，为要保证差分方程(4.9)的解是不增长的，的解是不增长的，只要选取只要选取h充分小，使充分小，使)10. 4(. 11 h在在 =h 的复平面上，这是以的复平面上，这是以(- -1,0)为圆心，为圆心，1为为半径的单位圆内部半径的单位圆内部. 称为欧拉法的称为欧拉法的

38、绝对稳定域绝对稳定域，一般，一般情形可由下面定义情形可由下面定义. 定义定义6 单步法单步法(4.1)用于解模型方程用于解模型方程y = y，若得，若得到的解到的解yn+1=E(h )yn，满足，满足|E(h )|1，则称方法，则称方法(4.1)是是绝对稳定绝对稳定的的. 在在 =h 的平面上的平面上, 使使|E(h )|1的变量的变量围成的区域，称为围成的区域，称为绝对稳定区域绝对稳定区域，它与实轴的交称，它与实轴的交称为为绝对稳定区间绝对稳定区间.上页上页上页上页上页上页下页下页下页下页下页下页对欧拉法对欧拉法E(h )=1+h ，其，其绝对稳定域绝对稳定域为为|1+h |1,绝对稳定区间

39、绝对稳定区间为为- -2h 0,在例在例5中中 =- -100,- -2- -100h0,即即0h2/100=0.02为为稳定区间稳定区间，在例，在例4中取中取h=0.025，故它是不稳定的，当取故它是不稳定的，当取h=0.005时它是稳定的时它是稳定的. 对二阶对二阶R- -K方法方法，解模型方程，解模型方程(4.1)可得到可得到,2)(121nnyhhy 故故.2)(1)(2 hhhE 绝对稳定域绝对稳定域由由|E(h )|1得到，于是可得得到，于是可得绝对稳定区间绝对稳定区间为为- -2h 0，即，即0h2/ .上页上页上页上页上页上页下页下页下页下页下页下页 类似可得三阶及四阶类似可得

40、三阶及四阶R- -K方法方法的的E(h )分别为分别为.!3)(!2)(1)(32 hhhhE .!4)(!3)(!2)(1)(432 hhhhhE 由由|E(h )|1可得到相应的可得到相应的绝对稳定域绝对稳定域. 当当 为实数时，为实数时，则得则得绝对稳定区间绝对稳定区间，它们分别为，它们分别为(图形复杂图形复杂p296) 三阶显式三阶显式R- -K方法方法： 四阶显式四阶显式R- -K方法方法：./51. 20, 051. 2 hh即即./78. 20, 078. 2 hh即即从以上讨论可知显式从以上讨论可知显式R- -K方法方法的的绝对稳定域绝对稳定域均为均为有限域，都对步长有限域，都

41、对步长h有限制有限制. 如果如果h不在所给的不在所给的绝对稳绝对稳定区间定区间内，方法就不稳定内，方法就不稳定.上页上页上页上页上页上页下页下页下页下页下页下页 例例5 分别取分别取h=0.1及及h=0.2，用经典的四阶，用经典的四阶R- -K方方法法(3.13)计算计算初值问题初值问题 . 1)0()10(,20yxyy 解解 本例本例 =- -20, ,h 分别为分别为- -2及及- -4，前者在，前者在绝对绝对稳定区间稳定区间内，后者则不在，用四阶内，后者则不在，用四阶R- -K方法方法计算其计算其误差见下表误差见下表xn0.20.40.60.81.0h=0.1h=0.20.9310-

42、-14.980.1210- -125.00.1410- -2125.00.1510- -3625.00.1710- -43125.0从以上结果看到，如果步长从以上结果看到，如果步长h不满足不满足绝对稳定条绝对稳定条件件，误差增长很快，误差增长很快. .上页上页上页上页上页上页下页下页下页下页下页下页对隐式单步法对隐式单步法, ,可以同样讨论方法的可以同样讨论方法的绝对稳定性绝对稳定性, ,例如对例如对后退欧拉法后退欧拉法，用它解模型方程可得，用它解模型方程可得,111nnyhy 故故.11)( hhE 由由|E(h )|1，这是以，这是以(1,0)为圆心，为圆心，1为半径的单位圆外部为半径的单

43、位圆外部. 故方法的故方法的绝对稳定区绝对稳定区间间为为- - h 0. 当当 0时，则时，则0h，即对任何步长均，即对任何步长均为稳定的为稳定的.上页上页上页上页上页上页下页下页下页下页下页下页 对对隐式梯形法隐式梯形法，它用于解模型方程，它用于解模型方程(4.8)得得,21211nnyhhy 故故.2121)( hhhE 对对Re( )0有有|E(h )|1，故，故绝对稳定域绝对稳定域为为 =h 的左半的左半平面，平面，绝对稳定区间绝对稳定区间为为- - h 0，即，即0h 时时隐式梯隐式梯形法形法均是稳定的均是稳定的.上页上页上页上页上页上页下页下页下页下页下页下页分析得到分析得到隐式欧

44、拉法隐式欧拉法与与梯形方法梯形方法的的绝对稳定域绝对稳定域均均为为h | Re(h )0，在具体计算中步长，在具体计算中步长h的选取只需考的选取只需考虑计算精度及迭代收敛性要求而不必考虑稳定性，具虑计算精度及迭代收敛性要求而不必考虑稳定性，具有这种特点的方法需特别重视有这种特点的方法需特别重视, 由此给出下面的定义由此给出下面的定义. 定义定义7 如果数值方法的绝对稳定域包含了集合如果数值方法的绝对稳定域包含了集合h | Re(h )0，那么称此方法是，那么称此方法是A- -稳定稳定的的. .由定义知由定义知A- -稳定稳定方法对步长方法对步长h没有限制没有限制. .上页上页上页上页上页上页下

45、页下页下页下页下页下页9.5 线性多步法线性多步法 在逐步推进的求解过程中在逐步推进的求解过程中, 计算计算yn+1之前事实上之前事实上已经求出了一系列的近似值已经求出了一系列的近似值y0,y1,yn, 如果充分利如果充分利用前面多步的信息来预测用前面多步的信息来预测yn+1, 则可以期望会获得较则可以期望会获得较高的精度高的精度. 这就是构造所得这就是构造所得线性多步法线性多步法的基本思想的基本思想. 构造构造多步法多步法的主要途径基于数值积分方法和基的主要途径基于数值积分方法和基于泰勒展开方法，前者可直接由方程于泰勒展开方法，前者可直接由方程(1.1)两端积分两端积分后利用插值求积公式得到

46、后利用插值求积公式得到. 本节主要介绍基于泰勒本节主要介绍基于泰勒展开的构造方法展开的构造方法.上页上页上页上页上页上页下页下页下页下页下页下页9.5.1 线性多步法的一般公式线性多步法的一般公式 如果计算如果计算yn+k时，除用时，除用yn+k- -1的值，还要用到的值，还要用到yn+i (i=0,1,k- -2)的值，则称此方法为的值，则称此方法为线性多步法线性多步法. 一般一般的的线性多步法公式线性多步法公式可表示为可表示为) 1 . 5(,010 kiinikiiniknfhyy 其中其中yn+1为为y(xn+1)的近似，的近似，fn+i=f(xn+i, yn+i), 这里这里xn+i

47、=xn+ih, i, i为常数为常数, 0及及 0不全为零不全为零, 则称则称(5.1)为为线性线性k步法步法, 计算时需先给出前面计算时需先给出前面k个近似值个近似值y0,y1,yk- -1, 再由再由(5.1)逐次求出逐次求出 yk,yk+1,.上页上页上页上页上页上页下页下页下页下页下页下页 如果如果 k=0，则，则(5.1)称为称为显式显式k步法步法，这时，这时yn+k可可直接由直接由(5.1)算出；如果算出；如果 k0, 则则(5.1)称为称为隐式隐式k步法步法,求解时与梯形法求解时与梯形法(2.7)相同相同, 要用迭代法方可算出要用迭代法方可算出yn+k. (5.1)中系数中系数

48、i及及 i可根据方法的局部截断误差及阶确可根据方法的局部截断误差及阶确定，其定义为定，其定义为 定义定义8 设设y(x)是初值问题是初值问题(1.1), (1.2)的准确解，的准确解，线性多步法线性多步法(5.1)在在xn+k上上局部截断误差局部截断误差为为);(hxyLTnkn ) 2 . 5(. )()()(010 kiinikiiniknxyhxyxy 若若Tn+k=O(hp+1)，则称方法，则称方法(5.1)是是p阶的阶的，如果，如果p 1, 则则称方法称方法(5.1)与微分方程与微分方程(1.1)是是相容的相容的.上页上页上页上页上页上页下页下页下页下页下页下页由定义由定义8，对，对

49、Tn+k在在xn处泰勒展开，由于处泰勒展开，由于,)(! 3)()(! 2)()()()(32 nnnnnxyihxyihxyihxyihxy.)(! 2)()()()(2 nnnnxyihxyihxyihxy代入代入(5.2)得得) 3 . 5(,)()()()()(2210 npppnnnknxyhcxyhcxyhcxycT上页上页上页上页上页上页下页下页下页下页下页下页其中其中., 3 , 2)2()!1(1) 1(2(!11110121 qkqkkqckqqkqqqq ),(11100 kc ),() 1(2101211kkkkc ) 4 . 5( 上页上页上页上页上页上页下页下页下页

50、下页下页下页若在公式若在公式(5.1)中选择系数中选择系数 i 及及 i ，使它满足，使它满足. 0, 0110 ppcccc由定义可知此时所构造的多步法是由定义可知此时所构造的多步法是p阶的，且阶的，且) 5 . 5().()(2) 1(11 pnpppknhOxyhcT称右端第一项为称右端第一项为局部截断误差主项局部截断误差主项, cp+1称为称为误差常数误差常数.根据相容性定义根据相容性定义, p 1，即，即c0=c1=0，由，由(5.4)得得) 6 . 5(., 1011110 kikiikiik 故线性多步法故线性多步法(5.1)与微分方程与微分方程(1.1)相容的充分必要条相容的充

51、分必要条件是件是(5.6)成立成立.上页上页上页上页上页上页下页下页下页下页下页下页显然，当显然，当k=1时，若时，若 1=0，则由，则由(5.6)可求得可求得 0=1， 0=1.此时公式此时公式(5.1)为为,1nnnhfyy 即为欧拉法即为欧拉法. 从从(5.4)可求得可求得c2=1/20，故方法为，故方法为1阶精阶精度，且局部截断误差为度，且局部截断误差为).()(21321hOxyhTnn 这和第这和第2节给出的定义及结果是一致的节给出的定义及结果是一致的.上页上页上页上页上页上页下页下页下页下页下页下页对对k=1，若，若 10，此时方法为隐式公式，为了确，此时方法为隐式公式，为了确定

52、系数定系数 0, 0, 1, 可由可由c0=c1=c2=0解得解得 0=1, 0= 1=1/2. 于是得到公式于是得到公式).(211 nnnnffhyy即为梯形公式即为梯形公式.由由(5.4)可求得可求得c2=- -1/12，故，故p=2，所以梯形法是二，所以梯形法是二阶方法，其局部截断误差主项是阶方法，其局部截断误差主项是).(123nxyh 这与这与9.2节中讨论是一致的节中讨论是一致的.对对k 2的多步法公式都可利用的多步法公式都可利用(5.4)确定系数确定系数 i, i,并由并由(5.5)式给出局部截断误差，下面只就若干常用的式给出局部截断误差，下面只就若干常用的多步法导出具体公式多

53、步法导出具体公式.上页上页上页上页上页上页下页下页下页下页下页下页9.5.2 阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式p299自学自学.9.5.3 米尔尼方法与辛普森方法米尔尼方法与辛普森方法p301自学自学.9.5.4 汉明方法汉明方法p302自学自学.9.5.5 预测预测-校正方法校正方法p303自学自学.9.5.6 构造多步法公式的注记和例构造多步法公式的注记和例p305自学自学.上页上页上页上页上页上页下页下页下页下页下页下页9.6 线性多步法的收敛性与稳定性线性多步法的收敛性与稳定性p306自看自看.上页上页上页上页上页上页下页下页下页下页下页下页9.7 一阶方程组与刚性方程组一

54、阶方程组与刚性方程组9.7.1 一阶方程组一阶方程组 前面我们研究了单个方程前面我们研究了单个方程y = =f 的数值解，只要的数值解，只要把把y 和和f 理解为向量，那么，所提供的各种计算公式理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形即可应用到一阶方程组的情形. 考察考察一阶方程组一阶方程组), 2 , 1(),(21NiyyyxfyNii 的初值问题，初始条件给为的初值问题，初始条件给为)., 2 , 1(,)(00Niyxyii 上页上页上页上页上页上页下页下页下页下页下页下页若采用向量的记号，记若采用向量的记号，记(向量向量),),(21TNyyyy ,),(00

55、2010TNyyyy .),(21TNffff 则上述方程组的初值问题可表示为则上述方程组的初值问题可表示为) 1 . 7(.)(),(00 yxyyxfy上页上页上页上页上页上页下页下页下页下页下页下页求解这一初值问题的求解这一初值问题的四阶龙格四阶龙格- -库塔公式库塔公式为为(向量向量),22(643211kkkkhyynn ).,(),2,2(),2,2(),(3423121hkyhxfkkhyhxfkkhyhxfkyxfknnnnnnnn式中式中(向量向量)上页上页上页上页上页上页下页下页下页下页下页下页或表示为或表示为(分量分量), 2 , 1()22(643211,NiKKKKhyyiiiiinni 其中其中(分量分量) ).,(),2,2,2,2(),2,2,2,2(),(323213142222121312121112211NNnnnniiNNnnnniiNNnnnniiNnnnniihKyhKyhKyhxfKKhyKhyKhyhxfKKhyKhyKhyhxfKyyyxfK这里这里yin是第是第i个因变量个因变量yi(x)在节点在节点xn=x0+nh的近似值的近似值.上页上页上页上页上页上页下页下页下页下页下页下页为了帮助理解这一公式的计算过程，