数值分析09-常微方程数值解法.ppt

《数值分析09-常微方程数值解法.ppt》由会员分享，可在线阅读，更多相关《数值分析09-常微方程数值解法.ppt（81页珍藏版）》请在三一办公上搜索。

1、阜师院数科院第八章 常微分方程数值解法,8-1,第八章,常微分方程数值解法,阜师院数科院第八章 常微分方程数值解法,8-2,第八章目录,1 欧拉（Euler）方法 1.1 Euler法及其简单改进 1.2 改进的Euler法2 龙格库塔（Runge-kutta）方法 2.1 龙格-库塔方法的基本思想 2.2 二阶龙格-库塔公式 2.3 高阶R-K公式 2.4 变步长R-K法3 线性多步法4 一阶方程组与高阶方程初值问题5 收敛性与稳定性,阜师院数科院第八章 常微分方程数值解法,8-3,第八章 序,许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动

2、，工程问题中的电路分析等，都可归结为常微分方程的初值问题。,所谓初值问题，是函数及其必要的导数在积分的起始点为已知的一类问题，一般形式为：,我们先介绍简单的一阶问题：,阜师院数科院第八章 常微分方程数值解法,8-4,第八章 序,由常微分方程的理论可知：上述问题的解唯一存在。常微分方程求解求什么？应求一满足初值问题（81）的解函数y=y(x)，如对下列微分方程：,高等数学中，微分方程求解，如对一阶微分方程：y=f(x,y)是求解解函数y=y(x)，使满足上述方程。但能够求出准确的解析函数y(x)的微分方程是很少的，高数中研究微分方程的求解，是分门别类讨论，对不同类型的微分方程，求解方法不一样，因

3、此，要求解微分方程，首先必须认清类型。,阜师院数科院第八章 常微分方程数值解法,8-5,微分方程 数值解,而常微分方程 初值问题的数值解法，是要寻求解函数y(x)在一系列点y(xi)（离散点）:,上 y(xi)的近似值yi（i=1,2,n），并且还可由这些（xi,yi)（i=1,2,n）构造插值函数作为近似函数。上述离散点相邻两点间的距离hi=xi-1-xi 称为步长，若hi 都相等为一定数h,则称为定步长，否则为变步长。,由于在实际问题和科学研究中遇到的微分方程往往很复杂，绝大多数很难，甚至不可能求出解析函数y(x)，因此只能考虑求其数值解。,本章重点讨论如下一阶微分方程：,在此基础上介绍一

4、阶微分方程组与高阶微分方程的数值解法。,阜师院数科院第八章 常微分方程数值解法,8-6,1 欧拉（Euler）法,以Euler法及其改进方法为例,说明常微分方程初值问题数值解法的一般概念，Euler法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其它方法的特点。,Euler法用于求解一阶微分方程初值问题：,阜师院数科院第八章 常微分方程数值解法,8-7,1.1 Euler法及其简单改进,Euler公式为：,由x0出发x1,x2,xN，而利用此式可算出对应的 y1,y2,yN，式（8-2）称为差分方程（序列yn满足的方

5、程）。,下面是Euler公式的推导：,一、从几何意义出发：y=f(x,y)的解函数y=y(x)在xoy平面上是一族解曲线，而初值问题则是其中一条积分曲线，假定y=y(x)的曲线如图8-1从给定的点P0(x0,y0)出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率 y=f(x0,y0)，因此可 得切线：（点斜式）,阜师院数科院第八章 常微分方程数值解法,8-8,Euler公式的推导(续1）,几何意义：用折线近似解曲线y=y(x)，折线不会偏离太远，因为每项以f(x,y)（斜率）修正。,切线与x=x1交于P1(x1,y1)，在x0,x1上以切线,近似曲线，,阜师院数科院第八章 常微分方

6、程数值解法,8-9,Euler公式的推导（续2）,二、利用Taylor级数：将y(x)在xn处展开：,阜师院数科院第八章 常微分方程数值解法,8-10,Euler公式的推导（续3）,公式（8-3）称为后退Euler公式,所谓局部载断误差是指以yn代替y(xn)而得到y(xn+1)的近似值yn+1的误差（只估计这一步的误差）。,三、利用数值微分公式：利用两点公式,后退Euler公式的局部截断误差为:,阜师院数科院第八章 常微分方程数值解法,8-11,Euler公式的推导（续4）,阜师院数科院第八章 常微分方程数值解法,8-12,Euler公式的推导（续5）,四、利用数值积分公式：在xn,xn+1

7、上对y(x)=f(x,y(x)积分,对右端积分项采用不同的数值积分公式，便可得到各种不同的求解dE初值问题的计算公式。如，以矩形面积代替曲边梯形面积,1）以左矩形面积代替曲边梯形面积如图8-2，亦即以,图8-2,阜师院数科院第八章 常微分方程数值解法,8-13,3）以梯形公式（面积）代替曲边形如图8-4则有,式（8-5）称为求dE初值问题的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数得到的，而Euler公式则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的解微分方程的数值公式，梯形公式比Euler公式更

8、准确一些，误差更小。,Euler公式的推导（续6）,2）以右矩形面积代替曲边梯形：如图8-3亦即以,阜师院数科院第八章 常微分方程数值解法,8-14,Euler公式注释,注1：Euler公式为显式，右矩形，梯形公式为隐式；,注2：Euler公式，梯形，右矩形公式为单步法，计算yn+1只用yn，而中点法公式为多步法（还可如上二所述，构造多步法）即必须已知yn-1，yn才 能计算yn+1，（求y0,y1不能用此公式。y0,y1称为多步法的开始值，y0给定，而y1必须由其它公式算出，然后才能用中点法）；,注3：前面已有Euler法 的局部截断误差：,后退Euler法的局部截断误差：,误差阶：如果局部

9、截断误差,则称方法为P阶的。,阜师院数科院第八章 常微分方程数值解法,8-15,显然，步长h越小，阶数P越高，局部截断误差越小，当然计算精度越高；,注4：梯形法是几阶？梯形法精度比Euler法高，阶数肯定比Euler法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积 分公式计算,因此由积分中梯形公式的误差知此,时的局部截断误差为：,梯形法为2阶方法！,Euler法，后退Euler法为1阶方法，而中点法为2 阶，,Euler公式注释（续）,阜师院数科院第八章 常微分方程数值解法,8-16,关于Euler法的整体截断误差注释,注5：关于Euler法的整体截断误差：,Euler方

10、法的局部截断误差公式为：,实际计算时，yn是y(xn)的近似值，因此，计算过程中除每步所产生的局部截断误差外，还有因前面的计算不准确而引起的误差。在不考虑舍入误差的情况下，称y(xn+1)与yn+1之差为整体截断误差，记为：,下面讨论Euler方法的整体截断误差。,为简便起见，假定函数f(x,y)充分光滑，问题（8-1）解y(x)在a,b上二阶连续可微，于是由式（8-6），局部截断误差有界，即存在M0，使得对任意xa,b，都有|y(x)|M，从而有：,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-17,式（8-7）表明，Euler方法的整体截断误差与h同阶，当h0时，en0。,关于

11、Euler法的整体截断误差注释（续）,反复递推得：,阜师院数科院第八章 常微分方程数值解法,8-18,结 论,对于实际问题来说，由于L，M 难以估计，因此（8-6）很难应用，而且上述推导过程中一再放大了误差上限，这样的估计往往也很保守，远远大于实际的误差，但是，从估计式（8-7）却可以得到下面很有用处的结论。,1）当h0时，en0即，,亦即数值解yn，一致收敛于初值问题（8-1）的真解y(xn)，,并且，Euler法的整体截断误差的阶为O(h)与h同阶，比局部截断误差低一阶。,2）舍入误差 局部截断误差 对实际计算结果有影响，并且随h减少 而减少或增大。,阜师院数科院第八章 常微分方程数值解法

12、,8-19,3）计算结果与解法的阶数p，真解的导数y(p+1)有关，p越大，h p+1越小，|y(p+1)()|的上限越大，M 也越大，因此为保证精度当然应选阶数p较高的方法。但如果M 很大，当f(x,y)是分段连续的函数时，则应采用低阶的方法如用Euler法。,结 论（续）,4）计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部 截断误差，对多步法，应采用跟多步法同阶的方法计算开始值。,阜师院数科院第八章 常微分方程数值解法,8-20,1.2 改进的Euler法,梯形公式为二阶方法，但却是隐式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，计算量较大，通常在实际计算时

13、，将Euler法与梯形公式合起来使用，即先使用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为预测值，然后用梯形公式去提高精度，将yn+1(0)校正为较准确的值：,由于函数f(x,y)满足Lipschitz条件，容易得出:,阜师院数科院第八章 常微分方程数值解法,8-21,改进的Euler法(续）,阜师院数科院第八章 常微分方程数值解法,8-22,预测校正型公式,实际经验表明，式（8-8）的迭代效果主要体现在第一次，由此构成如下的预测校正型公式:,此式称为改进的Euler公式，为上机计算编程方便，常将式（8-9）改写为:,下面分析改进的Euler公式的局部截断误差：,阜师

14、院数科院第八章 常微分方程数值解法,8-23,改进的Euler公式的局部截断误差分析,假定yn=y(xn)，y(xn+1)的Taylor展式为:,对于改进的Euler公式，由于,这说明改进的Euler法的局部截断误差为O(h3)，比Euler公式高一阶，是二阶方法。,阜师院数科院第八章 常微分方程数值解法,8-24,改进的Euler公式举例,例1,这些结果在表8-1中，可见计算结果的精度，Euler法与后退Euler法差不多，与准确值相比较Euler法偏小，而后退Euler法偏大；中点法与梯形法精度同为2阶，但梯形法更好一些，这跟它们局部截断误差的符号，阶数和系数的大小是完全一致的。,表见下屏

15、：,阜师院数科院第八章 常微分方程数值解法,8-25,表格8-1,表8-1 y=y,y(0)=1 的数值解(h=0.1),阜师院数科院第八章 常微分方程数值解法,8-26,表格8-2,而表8-2是分别取了不同的h=0.1,h=0.01，h=0.001,h=0.0001，还是利用这些公式，经过若干步的计算（h越小，计算量越大）算到y(1)的近似值，可见：随着h的减小，y(1)的近似值的精度在提高，0.01比0.001差，即0.001比0.01时的y(1)准确。,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-27,表8-2计算结果说明（续）,但h太小，到h=0.0001时却又变得误差大

16、了，这与前面所说h越小，p阶越高，应该局部截断误差越小，因而计算精度更高矛盾了，为什么会产生这种情况呢？这是由于h太小而引起计算量大因而造成了舍入误差和截断误差的积累，这种情况由于初值问题不同可能会影响更大，偏离更严重，如下面的例2。这种问题实际上是稳定性问题，我们将会讨论方法的稳定性，由此得出对h有一定的要求的稳定性制区域。,阜师院数科院第八章 常微分方程数值解法,8-28,例2,求解初值问题y=-20y，y(0)=1，计算y(1)的近似值。,解：类似于例1，用欧拉法、后退欧拉法、中点法、梯形法求解，得到如下表8.3。表8.3 y=20,y(0)=1的解y(1)的近似值(y(1)=0.206

17、116E8),由表8.3可见，尽管中点法的阶数与梯形法相同，比欧拉法和后退欧拉法的阶数高，计算结果的精度却很糟糕。此外，尽管欧拉法与后退欧拉法的阶数相同，但欧拉法计算结果的精度，当h=0.1时却比后退欧拉法差。,阜师院数科院第八章 常微分方程数值解法,8-29,2 龙格-库塔（Runge-kutta）方法,2.1 龙格-库塔方法的基本思想：,因此只要对平均斜率k*提供一种算法，由（8-11）式便相应地得到一种微分方程的数值计算公式。,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-30,龙格-库塔方法的基本思想(续）,改进欧拉公式比欧拉公式精度高的原因，也就在于确定平均斜率时，多取了

18、一个点的斜率值。因此它启发我们，如果设法在xi,xi+1上多预报几个点的斜率值，然后将它们加权平均作为k*的近似值，则有可能构造出更高精度的计算公式，这是龙格-库塔方法的基本思想。,用这个观点来研究欧拉公式与改进欧拉公式，可以发现欧拉公式由于仅取xn一个点的斜率值f(xn,yn)作为平均斜率k*的近似值，因此精度很低。而改进欧拉公式（8-10）却是利用了xn与xn-1两个点的斜率值k1=f(xn,yn)与k2=f(xn+1,yn+hk1)取算术平均作为平均斜率k*的近似值。,其中k2是通过已知信息yn利用欧拉公式求得的。,阜师院数科院第八章 常微分方程数值解法,8-31,2.2 二阶龙格-库塔

19、公式,首先推广改进欧拉公式，考察区间xn,xn+1内任一点：,我们希望用xn和x n+1两个点的斜率值k1和k2加权平均作为 平均斜率k*的近似值：,其中c1，c2为待定常数，同改进欧拉公式一样，这里仍取：,问题在于怎样预测xn+l处的斜率值k2。仿照改进欧拉公式，先用 欧拉公式提供y(xn+l)的预测值,然后再用预测值yn+l通过计算f 产生斜率值k2=f(xn+l,yn+l)，这样设计出的计算公式具有形式：,阜师院数科院第八章 常微分方程数值解法,8-32,二阶龙格-库塔公式（续1）,公式（8-12）中含有三个待定参数c1,c2和l，我们希望适当选取这些参数值，使得公式（8-12）具有二阶

20、精度，亦即使：,现在仍假定yn=y(xn)，即yn是准确的，将y(xn+1)与yn+1都在xi处作泰勒展开：,阜师院数科院第八章 常微分方程数值解法,8-33,二阶龙格-库塔公式（续2）,代入（8-12）式，得：,比较（8-13）与（8-14）两式，要使公式具有二阶精度，只有满足下列条件：,这里一共有三个待定参数，但只需满足两个条件，因此有一个自由度，于是满足条件（8-15）的参数不止一组，而是一族，相应的公式（8-12）也有一族，这些公式统称为二阶龙格-库塔公式，简称二阶R-K公式。,阜师院数科院第八章 常微分方程数值解法,8-34,特别，当l=1即xn+l=xn+1时，c1=c2=1/2，

21、二阶R-K公式就 是改进欧拉公式。如果取l=1/2，则c1=0,c2=1，这时二阶R-K公式称为变形的欧拉公式，其形式见左边：,从表面上看，变形的欧拉公式仅含一个斜率值k2，但k2是通过k1计算出来的，因此每完成一步，仍然需要两次计算函数f 的值，工作量和改进欧拉公式相同。,二阶龙格-库塔公式（续3）,阜师院数科院第八章 常微分方程数值解法,8-35,构造二阶R-K公式的主要步骤,综上所述，构造二阶R-K公式主要由以下几步产生：,在区间xn,xn+1上取二点，预报相应点 的斜率值；2)对此两斜率值加权平均作为平均斜率 值的近似值；3)将yn+1与y(xn+1)都在xi处作泰勒展开，为使公式达到

22、二阶精度，比较相应系 数，建立有关参数所应满足的方程组；4)解此方程组得一族二阶R-K公式。,阜师院数科院第八章 常微分方程数值解法,8-36,2.3 高阶R-K公式,为了进一步提高精度，在xn,xn+1上除xn和xn+l外再增加一点xn+m=xn+mh(l m 1)，并用xn,xn+l,xn+m三处的斜率值k1,k2,k3加权平均作为k*的近似值，这时计算公式为：,其中k1,k2仍用（8-12）式所取的形式。为了预测xn+m处的斜率值k3，要定出xn+m处所对应的yn+m，可以看作在区间xn,xn+m上使用二阶R-K公式，从而得到y(xn+m)的预测值：,于是，再通过计算函数值f 得到:,阜

23、师院数科院第八章 常微分方程数值解法,8-37,高阶R-K公式（续1）,这样设计出的计算公式具有形式:,运用泰勒展开方法，适当选择参数c1,c2,c3,l,m,及b1,b2使上述公式具有三阶精度采用与(8-12)类似的处理方法，得到这些参数需要满足的条件：,(见下屏),阜师院数科院第八章 常微分方程数值解法,8-38,高阶R-K公式（续2）,满足5个条件的一族公式（8-18），统称为三阶R-K公式，其中常见的是库塔公式：c1=1/6,c2=4/6,c3=1/6,l=1/2,m=1,b1=-1,b2=2,阜师院数科院第八章 常微分方程数值解法,8-39,若需再将精度提高至四阶，用上述类似处理方法

24、，只是必须在xi,xi+1上用四个点处的斜率加权平均作为k*的近似值，构成一族四阶龙格-库塔公式，由于推导复杂，这里从略，只将常用的公式介绍如下：,四阶公式中常用的还有下面的Gill公式：,高阶R-K公式（续3）,一般称它为标准四阶R-K公式，其局部截断误差是O(h5)。,阜师院数科院第八章 常微分方程数值解法,8-40,Gill 公式,在计算机上它具有节省存贮单元和控制舍入误差增长的优点。,阜师院数科院第八章 常微分方程数值解法,8-41,初值问题两种方法举例,用改进的Euler法和四阶标准R-K公式解初值问题：,例3,解 用改进的Euler公式，取h=0.1计算公式为：,部分计算结果见表8

25、-4：,表8-4,阜师院数科院第八章 常微分方程数值解法,8-42,对四阶标准R-K法，取h=0.2，计算公式如下：,例3（续1）,阜师院数科院第八章 常微分方程数值解法,8-43,例 3（续2）,表8-5,计算结果见表8-5，与准确解比较，至少有四位有效数字。而取步长h=0.1时，用改进的欧拉公式计算，最多只有三位有效数字（见表8-4），虽然改进的欧拉公式每前进一步只要计算两次f 值，而四阶R-K法要计算4次f 值，但由于后者步长比前者大，所耗费的计算总量还是差不多的。这个例子又一次显示了选择算法的重要意义。,阜师院数科院第八章 常微分方程数值解法,8-44,表8-6计算f 的次数与精度阶数

26、的关系,表8-6 的说明,一般，高精度的方法要求解有较好的光滑性，解的光滑性是由f(x,y)的可微性决定。如果解的光滑性差，则用四阶R-K法求得的数值解反而不如用改进的欧拉法好。因此，一定要针对具体问题的特点来选择求解方法。,例 3（续3）,阜师院数科院第八章 常微分方程数值解法,8-45,表8-6 的说明（续）,龙格-库塔方法的推导是在用泰勒展开方法的基础上进行的，因而它要求所求的解具有较好的光滑性质。假若解的光滑性差，那么使用四阶R-K公式求得的数值解，其精度可能反而不如改进欧拉公式。在实际计算时，应当针对问题的具体特点，选择合适的算法。,从理论上讲，可以构造任意高阶的龙格-库塔公式，但只

27、要注意到精度的阶数与计算函数值f(x,y)的次数之间的关系不是等量增加的（见表8-6），就可得出再提高公式阶数已没有多大意思了。因此更说明四阶R-K公式是兼顾了精度及计算量的较理想的计算公式。,需要指出的是：,阜师院数科院第八章 常微分方程数值解法,8-46,2.4 变步长R-K法,当用数值方法解微分方程时，单从一步来看，步长越小，局部截断误差就越小。但从整体来看，步长越小，在求解范围内所要完成的步数就会越多。这一方面增加了计算量；另一方面也增大了舍入误差的累积，因此有个如何合理选择步长的问题。在解决这个问题时，需要考虑两个问题：怎样衡量和检验计算结果的精度？如何依据获得的精度信息及时处理步长

28、。,以四阶标准R-K公式为例，从节点xn出发，先以步长h求得xn+1处解的近似值，记作yn+1h，由于公式的局部截断误差为O(h5)，故有：,阜师院数科院第八章 常微分方程数值解法,8-47,变步长R-K法（续1）,式中c与y(5)(x)在区间xn,xn+1内的值有关，假设y(5)(x)在区间xn,xn+1内变化不大，则可将c 视作常数。将步长折半，即以h/2为步长，跨两步到xn+1，再求得一个近似值,每跨一步的局部截断误差为c(h/2)5，因此,阜师院数科院第八章 常微分方程数值解法,8-48,具体做法,3.以h作基本步长，当 时，反复将步长折半，每折半一 次，由xn+1的步数增加一倍，直至

29、,2.如果，又不需要节点xn+h处的数值解值，则可将步 长加倍，只要保证所求节点处的y值符合精度要求即可。,阜师院数科院第八章 常微分方程数值解法,8-49,变步长四阶标准R-K法解初值问题举例,例4,解 取基本步长h=0.1，节点xk=kh，=106，按变步长四阶标准R-K法进行计算，其结果见表8-7，显然，尽管每一步的步长均只折半一次，但解的精度可达106，与准确解：,比较，说明以上分析是正确的。,表8-7见下屏,阜师院数科院第八章 常微分方程数值解法,8-50,表8-7,阜师院数科院第八章 常微分方程数值解法,8-51,3 线性多步法,前两节介绍的几种差分格式有一共同的特点，在计算yn+

30、1时仅用前一步求得的yn，而对前几步的信息未予启用，因而称之为单步法。现讨论多启用一些已知信息的线性多步法。线性多步法公式的一般形式为：,其中i、i均为与f,n无关的常数，r、r不同时为0。由于求yn+1时要到前r+1步的信息yn-r，yn及对应的f 值，因称它为r+1步法。-1=0，式（8-21）为显式，-10时，式（8-21）为隐式。,阜师院数科院第八章 常微分方程数值解法,8-52,3.1 阿当姆斯Adams公式,当式（8-21）中系数0=1,1=r=0时，称之为阿当姆斯（Adams）公式，这类公式很容易借助数值积分导出。设y(x)是初值问题（8-1）的准确解，则:,两边从xn到xn+1

31、积分，得:,对此式右端应用数值积分公式，则可导出不同的阿当姆斯公式。例如，记fj=f(xj,y(xj)，取数据点(xn-1,fn-1)，(xn,fn)构造f 的线性插值多项式:,其中：,阜师院数科院第八章 常微分方程数值解法,8-53,阿当姆斯Adams公式（续）,阜师院数科院第八章 常微分方程数值解法,8-54,r+1步阿当姆斯显式公式（又称Adams-Bashforth公式）,类似地，取数据点（xn-r,fn-r）、(xn-r+1,fn-r+1)、(xn,fn)，构造f(x,y(x)的r 次插值多项式，可导出r+1步阿当姆斯显式公式（又称Adams-Bashforth公式）及其局部截断误差

32、:,其系数见如下表8-8，显然具有r+1阶精度。,阜师院数科院第八章 常微分方程数值解法,8-55,r+1步阿当姆斯隐式公式（又称Adams-Moulton公式）,若改取数据点（xn-r+1,fn-r+1）、(xn-r+2,fnr+2)、(xn+1,fn+1)，构造f(x,y(x)的r次插值多项式，则可以导出n+1步阿当姆斯隐式公式（又称Adams-Moulton公式）及其局部截断误差：,阜师院数科院第八章 常微分方程数值解法,8-56,其系数见下表8-9，显然也具有r+1阶精度。,可见，一阶阿当姆斯显式公式（对应r=0）即为欧拉公式；一阶与二阶阿当姆斯隐式公式对应r=0与r=1，分别为后退欧

33、拉公式与梯形公式。,r+1步阿当姆斯隐式公式（又称Adams-Moulton公式）（续）,阜师院数科院第八章 常微分方程数值解法,8-57,四阶阿当姆斯显式公式,高阶阿当姆斯公式中最常用的是四阶公式，对应r=3。若选xn-3,xn-2，xn-1,xn作为插值节点，则四阶阿当姆斯显式公式（又称为外推公式）及其局部截断误差为：,此公式的优点是精度高，其局部截断误差与四阶R-K法同阶，但每前进一步只要计算一次函数值f(xn,yn)，其余的值f 在前面的计算中已求出，计算量远小于四阶R-K法。不足之处是不能从n=1开始使用，必须用其它方法（例如四阶R-K法）求得出发值y1、y2、y3。此外，若要中途改

34、变 步长也较麻烦。,阜师院数科院第八章 常微分方程数值解法,8-58,四阶阿当姆隐式公式内插公式及其局部截断误差为（取xn-2,xn-1,xn,xn+1为插值节点）：,它也具有公式（8-26）的优缺点，只是对yn+1必须提供一个预测值，然后通过迭代才能求得yn+1。将公式（8-26）与公式（8-27）联合使用，可组成如下预测一校正系统：,四阶阿当姆斯显式公式（续）,阜师院数科院第八章 常微分方程数值解法,8-59,阿当姆斯公式举例,例5,分别用四阶阿当姆斯显式公式（8-26）与阿当姆斯预测一校正公式（8-28）解例3中的初值问题。,解 取步长h=0.1，由于公式（8-9）、（8-28）均为四步

35、法，因此必须用其它方法求出发值。现采用四阶标准龙格-库塔法求y1,y2,y3 然后再分别用指定的方法进行计算，计算结果见表8-10。,表8-10,阜师院数科院第八章 常微分方程数值解法,8-60,阿当姆斯预测校正公式,为了进一步提高上述预测一校正法的计算精度，可以利用预测公式与校正公式同阶的特点，导出局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计值进行修正，以此来提高解的精度。由式（8-26）、（8-27）可知，公式（8-28）中预测公式与校正公式的局部截断误差分别为：,若记预测值为pn+1，校正值为cn+1，则：,紧接下屏,阜师院数科院第八章 常微分方程数值解法,8-61

36、,于是得pn+1与cn+1的事后误差估计：,阿当姆斯预测校正公式（续）,分别代替pn+1与cn+1能提高解的精度。为简单计，在未求出cn+1以前，用pncn代替pn+1cn+1。,阜师院数科院第八章 常微分方程数值解法,8-62,阿当姆斯带误差修正的预测校正公式,据此，得到下述带误差修正的预测一校正公式：,式中初始值p3与c3可取为0。,阜师院数科院第八章 常微分方程数值解法,8-63,3.2 哈明(Hamming)公式,一般的线性多步公式（8-21）还可用待系数法并借助泰勒展开导出。例如考虑确定如下形式的四阶公式：,其中yn表示fn。只需求系数j、j，使其局部截断误差T=O(h5)即可。局部

37、截断误差为：,将其在xn点进行泰勒展开，得：,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-64,哈明(Hamming)公式（续1）,令hk(k=0,1,4)的系数为0,得关于j、j的方程组:,方程组有六个未知数，而只有五个方程，因此有无穷多解，取1为自由参量，则有:,阜师院数科院第八章 常微分方程数值解法,8-65,哈明(Hamming)公式（续2）,当1=0时，得公式：,这是一个四阶隐式公式。,完全类似，可导出如下形式的四阶显式公式：,其中最常用的是：,阜师院数科院第八章 常微分方程数值解法,8-66,Hamming公式（预测一校正公式）,用显式公式（8-34）与隐式公式（8-

38、33）联合可组成预测一校正公式,此式称为哈明（Hamming）公式，大量数值实验表明，此公式的数值稳定性在同类公式中几乎是最好的。为了进一步改善哈明公式计算结果的精度而又不致增加过多的计算量，可仿照带误差修正的阿当姆斯预测一校正公式，建立一个带误差修正的哈明公式：,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-67,前面介绍的四套预测一校正公式（8-28）、（8-29）、（8-32）、（8-33）都是实际计算中常用的方法，它们的 共同特点是计算量省而精度高，但都必须用别的单步法 为其提供出发值，由于解题过程是递推进行的，所以出 发值的精度对以后的计算影响很大，故一般采用四阶龙 格一

39、库塔法求其出发值。,Hamming公式（预测一校正公式）（续）,阜师院数科院第八章 常微分方程数值解法,8-68,4 一阶方程组与高阶方程初值问题,4.1 一阶方程组,下面以两个方程的情形为例，介绍一阶方程组的数值解法。设初值问题：,若采用向量记号，记：,则初值问题（8-34）可表示为：,它与初值问题（8-1）形式上完全一致，前面介绍的解初值问题（8-1）的各种数值方法均可用于解一阶方程组初值问题（8-35），只要将公式中的y与f换成向量y与f即可。,阜师院数科院第八章 常微分方程数值解法,8-69,一阶方程组（续1）,例如：解初值问题（8-35）的欧拉公式为：,即：,写成分量形式即为：,又如

40、解初值问题（8-35）的四阶标准龙格-库塔公式为：,阜师院数科院第八章 常微分方程数值解法,8-70,一阶方程组（续2）,写成分量形式为：,阜师院数科院第八章 常微分方程数值解法,8-71,4.2 高阶方程初值问题,高阶方程初值问题常转化为一阶方程组求解。例如：,引进新变量z=y，则此问题等价于,对其使用四阶标龙格-库塔法，计算公式为,阜师院数科院第八章 常微分方程数值解法,8-72,高阶方程初值问题举例,例5,用四阶龙格-库塔方法在0，1上取步长h=0.2，求解二阶方程初值问题：,解 先将二阶方程化为一 阶方程组。令z=y，则得方程组：,使用四阶龙格-库塔公式，其相应的形式为：,其中,阜师院

41、数科院第八章 常微分方程数值解法,8-73,例5（续）,计算结果列于表8-11中：,与准确解y=e x比较，可知所得数值解除x=1外，均有五位有效数。,阜师院数科院第八章 常微分方程数值解法,8-74,5 收敛性与稳定性,通过前面的讨论可以看到，微分方程数值解法的基本思想是：通过某种离散化手段，将微分方程转化为差分格式求解。从理论上说，这里需要解决两个问题，一是当步长h0时，差分解，这就是差分格式的收敛性问题；二是利用差分格式求解时，初始误差及计算过程中的舍入误差能否得到控制，这就是差分格式的稳定性问题。本节对这两个问题作一简单介绍。,阜师院数科院第八章 常微分方程数值解法,8-75,5.1

42、收敛性,用某种差分格式解微分方程时，若对求解区间a,b中的任一点x，,当h0时差分解yn一致收敛于微分方程之解y(x)，,对x一致成立，则称该差分格式是收敛的。,显然，只有收敛的差分格式才有实用价值。可以证明，前几节介绍的单步法与多步法对满足解的存在唯一性定理条件的微分方程都是收敛的。,对于收敛的差分格式，从理论上说只要h充分小，差分解yn的精度就可以任意高。,阜师院数科院第八章 常微分方程数值解法,8-76,5.2 稳定性,实际计算中舍入误差几乎不可避免，这类误差在递推计算过程中会不会恶性增长，以致淹没了差分方程的真解，这是微分方程数值解法中的另一个重问题。由于研究的对象不同，关于稳定性的定

43、义很多，这里只介绍其中一种绝对稳定性。由于实际问题中微分方程（8-1）的右端f(x,y)复杂多样，这给研究稳定性带来了困难，因此通常取以下模型方程为基础进行研究。,其中为常数，0表示微分方程本身是稳定的。,若取步长h，用某种差分格式解模型方程时，任何一步产生的舍入误差在以后的计算中引起的误差逐渐减弱，则称该差分格式关于z=h是绝对稳定的，使差分格式稳定的所有z值的集合称为绝对稳定区间。,阜师院数科院第八章 常微分方程数值解法,8-77,稳定性（续1）,例如，将欧拉公式用于模型方程y=y，得,不妨设在yn上有误差n，则由此引起yn+1的误差n+1满足:,显然，要保证误差逐渐减弱，必须使z=h满足

44、:,阜师院数科院第八章 常微分方程数值解法,8-78,表8-12列出了几个常用公式的绝对稳定区间。从表中可以看到，隐式公式比同阶显式公式的稳定区间大，这就是为什么最好立足于隐式公式求解的原因。,稳定性（续1）,阜师院数科院第八章 常微分方程数值解法,8-79,稳定性（续3）,研究一般方程的稳定性时，相当于模型方程中的，例如对欧拉公式:,设yn有误差n，则由此引起yn+1的误差n+1满足:,所以有：,（紧接下屏）,阜师院数科院第八章 常微分方程数值解法,8-80,综上所述，当用数值方法解微分方程时，步长h的选择由两个因素确定。一是精度要求，由收敛性可知，h越小精度越高。二是稳定性要求，应落在所选公式的绝对稳定域中。虽然h十分小时两个条件都满足，但h太小计算步数势必大增，每步的微小误差累积起来也很可观。因此，实际计算中在满足精度要求与稳定性的条件下，步长h应尽可能取得大些。,式中n在yn+n与yn之间，与式（8-37）相比，可知 起着的作用，因此为了保证数值方法的绝对稳定，应落在它的绝对稳定域中。,稳定性（续4）,阜师院数科院第八章 常微分方程数值解法,8-81,第八章,结 束,