分析常微方程数值解法

分析常微方程数值解法第八章常微分方程数值解法8-1第一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-2第八章目录§1欧拉（Euler）方法

1.1Euler法及其简单改进1.2改进的Euler法§2龙格库塔（Runge-kutta）方法

2.1龙格-库塔方法的基本思想2.2二阶龙格-库塔公式2.3高阶R-K公式2.4变步长R-K法§3线性多步法§4一阶方程组与高阶方程初值问题§5收敛性与稳定性第二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-3第八章序

许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。

所谓初值问题，是函数及其必要的导数在积分的起始点为已知的一类问题，一般形式为：

我们先介绍简单的一阶问题：第三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-4第八章序由常微分方程的理论可知：上述问题的解唯一存在。

常微分方程求解求什么？应求一满足初值问题（8—1）的解函数y=y(x)，如对下列微分方程：《高等数学》中，微分方程求解，如对一阶微分方程：y=f(x,y)是求解解函数y=y(x)，使满足上述方程。但能够求出准确的解析函数y(x)的微分方程是很少的，《高数》中研究微分方程的求解，是分门别类讨论，对不同类型的微分方程，求解方法不一样，因此，要求解微分方程，首先必须认清类型。第四页，共八十一页，2022年，8月28日第八章常微分方程数值解法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)，因此只能考虑求其数值解。本章重点讨论如下一阶微分方程：在此基础上介绍一阶微分方程组与高阶微分方程的数值解法。第五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-6§1欧拉（Euler）法

以Euler法及其改进方法为例,说明常微分方程初值问题数值解法的一般概念，Euler法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其它方法的特点。Euler法用于求解一阶微分方程初值问题：第六页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-71.1Euler法及其简单改进

Euler公式为：由x0出发x1,x2,…,xN，而利用此式可算出对应的y1,y2,…,yN，式（8-2）称为差分方程（序列{yn}满足的方程）。下面是Euler公式的推导：一、从几何意义出发：y=f(x,y)的解函数y=y(x)在xoy平面上是一族解曲线，而初值问题则是其中一条积分曲线，假定y=y(x)的曲线如图8-1从给定的点P0(x0,y0)出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率

y

=f(x0,y0)，因此可得切线：（点斜式）P1P2y(x)P0

x2x1x0第七页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-8Euler公式的推导(续1）几何意义：用折线近似解曲线y=y(x)，折线不会偏离太远，因为每项以f(x,y)（斜率）修正。切线与x=x1交于P1(x1,y1)，在[x0,x1]上以切线近似曲线，

第八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-9Euler公式的推导（续2）二、利用Taylor级数：将y(x)在xn处展开：

第九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-10Euler公式的推导（续3）公式（8-3）称为后退Euler公式

所谓局部载断误差是指以yn代替y(xn)而得到y(xn+1)的近似值yn+1的误差（只估计这一步的误差）。三、利用数值微分公式：利用两点公式并且Euler公式的局部截断误差为:后退Euler公式的局部截断误差为:第十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-11Euler公式的推导（续4）第十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-12Euler公式的推导（续5）四、利用数值积分公式：在[xn,xn+1]上对y(x)=f(x,y(x))积分

对右端积分项采用不同的数值积分公式，便可得到各种不同的求解dE初值问题的计算公式。如，以矩形面积代替曲边梯形面积1）以左矩形面积代替曲边梯形面积如图8-2，亦即以yf(x,y)xnxxn+1图8-2

第十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-13yf(x,y)xnxxn+1图8-3yf(x,y)xnxxn+1图8-43）以梯形公式（面积）代替曲边形如图8-4则有

式（8-5）称为求dE初值问题的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数得到的，而Euler公式则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的解微分方程的数值公式，梯形公式比Euler公式更准确一些，误差更小。Euler公式的推导（续6）2）以右矩形面积代替曲边梯形：如图8-3亦即以

第十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-14Euler公式注释注1：Euler公式为显式，右矩形，梯形公式为隐式；注2：Euler公式，梯形，右矩形公式为单步法，计算yn+1只用yn，而中点法公式为多步法（还可如上二所述，构造多步法）即必须已知yn-1，yn才能计算yn+1，（求y0,y1不能用此公式。y0,y1称为多步法的开始值，y0给定，而y1必须由其它公式算出，然后才能用中点法）；注3：前面已有Euler法的局部截断误差：后退Euler法的局部截断误差：

误差阶：如果局部截断误差则称方法为P阶的。第十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-15显然，步长h越小，阶数P越高，局部截断误差越小，当然计算精度越高；

注4：梯形法是几阶？梯形法精度比Euler法高，阶数肯定比Euler法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积分公式计算因此由积分中梯形公式的误差知此时的局部截断误差为：∴梯形法为2阶方法！Euler法，后退Euler法为1阶方法，而中点法为2阶，Euler公式注释（续）第十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-16关于Euler法的整体截断误差注释注5：关于Euler法的整体截断误差：Euler方法的局部截断误差公式为：

实际计算时，yn是y(xn)的近似值，因此，计算过程中除每步所产生的局部截断误差外，还有因前面的计算不准确而引起的误差。在不考虑舍入误差的情况下，称y(xn+1)与yn+1之差为整体截断误差，记为：下面讨论Euler方法的整体截断误差。

为简便起见，假定函数f(x,y)充分光滑，问题（8-1）解y(x)在[a,b]上二阶连续可微，于是由式（8-6），局部截断误差有界，即存在M>0，使得对任意x[a,b]，都有|y′(x)|≤M，从而有：（紧接下屏）第十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-17式（8-7）表明，Euler方法的整体截断误差与h同阶，当h0时，en0。

关于Euler法的整体截断误差注释（续）

反复递推得：第十七页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-18结论

对于实际问题来说，由于L，M难以估计，因此（8-6）很难应用，而且上述推导过程中一再放大了误差上限，这样的估计往往也很保守，远远大于实际的误差，但是，从估计式（8-7）却可以得到下面很有用处的结论。

1）当h0时，en0即，亦即数值解yn，一致收敛于初值问题（8-1）的真解y(xn)，并且，Euler法的整体截断误差的阶为O(h)与h同阶，比局部截断误差低一阶。2）舍入误差局部截断误差对实际计算结果有影响，并且随h减少而减少或增大。第十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-193）计算结果与解法的阶数p，真解的导数y(p+1)有关，p越大，hp+1越小，|y(p+1)(ξ)|的上限越大，M也越大，因此为保证精度当然应选阶数p较高的方法。但如果M很大，当f(x,y)是分段连续的函数时，则应采用低阶的方法如用Euler法。结论（续）

4）计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部截断误差，对多步法，应采用跟多步法同阶的方法计算开始值。第十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-201.2改进的Euler法

梯形公式为二阶方法，但却是隐式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，计算量较大，通常在实际计算时，将Euler法与梯形公式合起来使用，即先使用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为

预测值，然后用梯形公式去提高精度，将yn+1(0)

校正为较准确的值：由于函数f(x,y)满足Lipschitz条件，容易得出:第二十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-21改进的Euler法(续）第二十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-22预测──校正型公式

实际经验表明，式（8-8）的迭代效果主要体现在第一次，由此构成如下的预测──校正型公式:此式称为改进的Euler公式，为上机计算编程方便，常将式

（8-9）改写为:下面分析改进的Euler公式的局部截断误差：

第二十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-23改进的Euler公式的局部截断误差分析假定yn=y(xn)，y(xn+1)的Taylor展式为:对于改进的Euler公式，由于

这说明改进的Euler法的局部截断误差为O(h3)，比Euler公式高一阶，是二阶方法。第二十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-24改进的Euler公式举例例1

这些结果在表8-1中，可见计算结果的精度，Euler法与后退Euler法差不多，与准确值相比较Euler法偏小，而后退Euler法偏大；中点法与梯形法精度同为2阶，但梯形法更好一些，这跟它们局部截断误差的符号，阶数和系

数的大小是完全一致的。表见下屏：第二十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-25表格8-1表8-1

y=y,y(0)=1的数值解(h=0.1)

x精确解欧拉法后退欧拉中点法梯形法.1.904837.900000.909091.900000.904762.2.808731.810000.826446.820000.818594.3.740818.729000.751315.736000.740633.4.670320.656100.683013.627800.670096.5.060531.590490.620921.601440.606278.6.548812.531441.654474.552512.548537.7.496585.478298.5131458.490938.496295.8.449329.430467.466507.454324.449029.9.406570.387421.424098.400073.4062641.367879.348679.385543.374310.367573第二十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法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)准确。Y′=-y,y(0)=1的解y(1)的近似值(y(1)=0.367879)h欧拉法后退欧拉法中点法梯形法0.1.348678.385543.374310.3675730.01.366033.369711.367944.3678770.001.367700.368052.367879.3678760.0001.367800.367800.367881.368020（紧接下屏）第二十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-27表8-2计算结果说明（续）

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

表8.3y=20,y(0)=1的解y(1)的近似值(y(1)=0.206116E8)h欧拉法后退欧拉法中点法梯形法.11.169351E-4.514229E+60.01.203704E-9.120746E-7.413244E+7.192743E-8.001.168300E-8.251090E-8.484136E+5.205979E-8.0001.202081E-8.210331E-8.475130+3.206103E-8

由表8.3可见，尽管中点法的阶数与梯形法相同，比欧拉法和后退欧拉法的阶数高，计算结果的精度却很糟糕。此外，尽管欧拉法与后退欧拉法的阶数相同，但欧拉法计算结果的精度，当h=0.1时却比后退欧拉法差。第二十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-29§2龙格-库塔（Runge-kutta）方法2.1龙格-库塔方法的基本思想：

因此只要对平均斜率k*提供一种算法，由（8-11）式便相应地得到一种微分方程的数值计算公式。（紧接下屏）第二十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-30龙格-库塔方法的基本思想(续）

改进欧拉公式比欧拉公式精度高的原因，也就在于确定平均斜率时，多取了一个点的斜率值。因此它启发我们，如果设法在[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利用欧拉公式求得的。第三十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-312.2二阶龙格-库塔公式

首先推广改进欧拉公式，考察区间[xn,xn+1]内任一点：我们希望用xn和xn+1两个点的斜率值k1和k2加权平均作为平均斜率k*的近似值：其中c1，c2为待定常数，同改进欧拉公式一样，这里仍取：问题在于怎样预测xn+l处的斜率值k2。仿照改进欧拉公式，先用欧拉公式提供y(xn+l)的预测值

然后再用预测值yn+l通过计算f产生斜率值k2=f(xn+l,yn+l)，这样设计出的计算公式具有形式：第三十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-32二阶龙格-库塔公式（续1）公式（8-12）中含有三个待定参数c1,c2和l，我们希望适当选取这些参数值，使得公式（8-12）具有二阶精度，亦即使：现在仍假定yn=y(xn)，即yn是准确的，将y(xn+1)与yn+1都在xi处作泰勒展开：第三十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-33二阶龙格-库塔公式（续2）代入（8-12）式，得：比较（8-13）与（8-14）两式，要使公式具有二阶精度，只有满足下列条件：

这里一共有三个待定参数，但只需满足两个条件，因此有一个自由度，于是满足条件（8-15）的参数不止一组，而是一族，相应的公式（8-12）也有一族，这些公式统称为二阶龙格-库塔公式，简称二阶R-K公式。第三十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-34

特别，当l=1即xn+l=xn+1时，c1=c2=1/2，二阶R-K公式就是改进欧拉公式。

如果取l=1/2，则c1=0,c2=1，这时二阶R-K公式称为变形的欧拉公式，其形式见左边：

从表面上看，变形的欧拉公式仅含一个斜率值k2，但k2是通过k1计算出来的，因此每完成一步，仍然需要两次计算函数f的值，工作量和改进欧拉公式相同。二阶龙格-库塔公式（续3）第三十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-35构造二阶R-K公式的主要步骤

综上所述，构造二阶R-K公式主要由以下几步产生：

在区间[xn,xn+1]上取二点，预报相应点

的斜率值；2)对此两斜率值加权平均作为平均斜率

值的近似值；3)将yn+1与y(xn+1)都在xi处作泰勒展开，

为使公式达到二阶精度，比较相应系

数，建立有关参数所应满足的方程组；4)解此方程组得一族二阶R-K公式。第三十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-362.3高阶R-K公式

为了进一步提高精度，在[xn,xn+1]上除xn和xn+l外再增加一点xn+m=xn+mh(l

m1)，并用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得到:第三十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-37高阶R-K公式（续1）这样设计出的计算公式具有形式:运用泰勒展开方法，适当选择参数c1,c2,c3,l,m,及b1,b2使上述公式具有三阶精度采用与(8-12)类似的处理方法，得到这些参数需要满足的条件：(见下屏)第三十七页，共八十一页，2022年，8月28日第八章常微分方程数值解法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第三十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-39

若需再将精度提高至四阶，用上述类似处理方法，只是必须在[xi,xi+1]上用四个点处的斜率加权平均作为k*的近似值，构成一族四阶龙格-库塔公式，由于推导复杂，这里从略，只将常用的公式介绍如下：四阶公式中常用的还有下面的Gill公式：高阶R-K公式（续3）一般称它为标准四阶R-K公式，其局部截断误差是O(h5)。第三十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-40Gill公式

在计算机上它具有节省存贮单元和控制舍入误差增长的优点。第四十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-41初值问题两种方法举例用改进的Euler

法和四阶标准

R-K公式解初值

问题：例3[解]用改进的Euler公式，取h=0.1计算公式为：部分计算结果见表8-4：

表8-4xn改进的欧拉法准确解y(xn)误差xn改进的欧拉法准确解y(xn)误差01.100.61.4859561.4832400.0027160.21.1840961.1832160.0000880.81.6164761.6124520.0040240.41.3433601.3416410.0017191.01.7378691.7320510.005818第四十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-42对四阶标准R-K法，取h=0.2，计算公式如下：例3（续1）第四十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-43例3（续2）表8-5

xn四阶R-K法误差xn四阶R-K法误差0110.61.4832810.0000410.21.1832290.0000130.81.6125130.0000610.41.3416670.000261.01.7321400.000089

计算结果见表8-5，与准确解比较，至少有四位有效数字。而取步长h=0.1时，用改进的欧拉公式计算，最多只有三位有效数字（见表8-4），虽然改进的欧拉公式每前进一步只要计算两次f值，而四阶R-K法要计算4次f值，但由于后者步长比前者大，所耗费的计算总量还是差不多的。这个例子又一次显示了选择算法的重要意义。第四十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-44表8-6

计算f的次数与精度阶数的关系

每步计算f的次数23456789精度的阶数23445667表8-6的说明

一般，高精度的方法要求解有较好的光滑性，解的光滑性是由f(x,y)的可微性决定。如果解的光滑性差，则用四阶R-K法求得的数值解反而不如用改进的欧拉法好。因此，一定要针对具体问题的特点来选择求解方法。例3（续3）第四十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-45表8-6的说明（续）

龙格-库塔方法的推导是在用泰勒展开方法的基础上进行的，因而它要求所求的解具有较好的光滑性质。假若解的光滑性差，那么使用四阶R-K公式求得的数值解，其精度可能反而不如改进欧拉公式。在实际计算时，应当针对问题的具体特点，选择合适的算法。从理论上讲，可以构造任意高阶的龙格-库塔公式，但只要注意到精度的阶数与计算函数值f(x,y)的次数之间的关系不是等量增加的（见表8-6），就可得出再提高公式阶数已没有多大意思了。因此更说明四阶R-K公式是兼顾了精度及计算量的较理想的计算公式。

需要指出的是：第四十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-462.4变步长R-K法

当用数值方法解微分方程时，单从一步来看，步长越小，局部截断误差就越小。但从整体来看，步长越小，在求解范围内所要完成的步数就会越多。这一方面增加了计算量；另一方面也增大了舍入误差的累积，因此有个如何合理选择步长的问题。在解决这个问题时，需要考虑两个问题：

①怎样衡量和检验计算结果的精度？②如何依据获得的精度信息及时处理步长。

以四阶标准R-K公式为例，从节点xn出发，先以步长h求得xn+1处解的近似值，记作yn+1[h]，由于公式的局部截断误差为O(h5)，故有：第四十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法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，因此

第四十七页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-48具体做法

1.若<

且相差不多，则取，并以h为步长计算xn+2对应的yn+2。3.以h作基本步长，当>

时，反复将步长折半，每折半一

次，由xn+1的步数增加一倍，直至2.如果<<

，又不需要节点xn+h处的数值解值，则可将步

长加倍，只要保证所求节点处的y值符合精度要求即可。第四十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-49变步长四阶标准R-K法解初值问题举例

例4解取基本步长h=0.1，节点xk=kh，=106，按变步长四阶标准R-K法进行计算，

其结果见表8-7，显然，尽管每一步的步长均只折半一次，但解的精度可达106，与准确解：比较，说明以上分析是正确的。表8-7见下屏第四十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-50步长折半次数xkyk准确解y(xk)10.10.1049920.1049917110.20.2198680.2198680710.30.3443410.3443407710.40.4779530.4779534810.50.6201150.620114510.60.7701350.7701352810.70.9272700.9272702410.81.0907541.090753610.91.2598311.259830411.01.4337811.4337808表8-7第五十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-51§3线性多步法

前两节介绍的几种差分格式有一共同的

特点，在计算yn+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）为隐式。第五十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-523.1阿当姆斯Adams公式

当式（8-21）中系数0=1,

1=…=r=0时，称之为阿当姆斯（Adams）公式，这类公式很容易借助数值积分导出。设y(x)是初值问题（8-1）的准确解，则:两边从xn到xn+1积分，得:对此式右端应用数值积分公式，则可导出不同的阿当姆斯公式。例如，记fj=f(xj,y(xj))，取数据点(xn-1,fn-1)，(xn,fn)构造f的线性插值多项式:其中：第五十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-53阿当姆斯Adams公式（续）第五十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-54r+1步阿当姆斯显式公式

（又称Adams-Bashforth公式）

类似地，取数据点（xn-r,fn-r）、(xn-r+1,fn-r+1)、…、(xn,fn)，构造f(x,y(x))的r次插值多项式，可导出r+1步阿当姆斯显式公式（又称Adams-Bashforth公式）及其局部截断误差:rArBr0Br1Br2Br3Br4r+1011

1/21231

5/1221223165

3/83245559379

251/7204720190127742616127425195/288其系数见如下表8-8，显然具有r+1阶精度。第五十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-55r+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公式）及其局部截断误差：第五十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-56

其系数见下表8-9，显然也具有r+1阶精度。

rArBr-1*Br0*Br1*Br2*Br3*r+1*011

1/21211

1/12212581

1/2432491951

19/7204720251646264106193/160可见，一阶阿当姆斯显式公式（对应r=0）即为欧拉公式；一阶与二阶阿当姆斯隐式公式对应r=0与r=1，分别为后退欧拉公式与梯形公式。r+1步阿当姆斯隐式公式

（又称Adams-Moulton公式）（续）

第五十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法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。此外，若要中途改变步长也较麻烦。第五十七页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-58四阶阿当姆隐式公式内插公式及其局部截断误差为（取xn-2,xn-1,xn,xn+1为插值节点）：

它也具有公式（8-26）的优缺点，只是对yn+1必须提供一个预测值，然后通过迭代才能求得yn+1。将公式（8-26）与公式（8-27）联合使用，可组成如下预测一校正系统：

四阶阿当姆斯显式公式（续）第五十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-59阿当姆斯公式举例例5分别用四阶阿当姆斯显式公式（8-26）与阿当姆斯预测一校正公式（8-28）解例3中的初值问题。解

取步长h=0.1，由于公式（8-9）、（8-28）均为四步法，因此必须用其它方法求出发值。现采用四阶标准龙格-库塔法求y1,y2,y3‘然后再分别用指定的方法进行计算，计算结果见表8-10。表8-10xnyn四阶标准R-K法显式公式（8-25）预测-校正公式（8-28）准确解0.11.095446

1.0954451150.21.183217

1.1832159560.31.264912

1.2649110460.4

1.3415521.3416411.3416407860.5

1.4140461.4142131.4142135620.6

1.4830191.4832391.4832396870.7

1.5489191.5491921.5481833380.8

1.6121161.6124501.6124515490.9

1.6729171.6733781.633200531.0

1.7315691.7320481.732050807第五十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-60阿当姆斯预测—校正公式

为了进一步提高上述预测一校正法的计算精度，可以利用预测公式与校正公式同阶的特点，导出局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计值进行修正，以此来提高解的精度。由式（8-26）、（8-27）可知，公式（8-28）中预测公式与校正公式的局部截断误差分别为：若记预测值为pn+1，校正值为cn+1，则：紧接下屏第六十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-61于是得pn+1与cn+1的事后误差估计：

阿当姆斯预测—校正公式（续）分别代替pn+1与cn+1能提高解的精度。为简单计，在未求出cn+1以前，用pncn代替pn+1cn+1。第六十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-62阿当姆斯带误差修正的预测—校正公式据此，得到下述带误差修正的预测一校正公式：式中初始值p3与c3可取为0。

第六十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-633.2哈明(Hamming)公式

一般的线性多步公式（8-21）还可用待系数法并借助泰勒展开导出。例如考虑确定如下形式的四阶公式：其中yn表示fn。只需求系数j、j，使其局部截断误差T=O(h5)即可。局部截断误差为：将其在xn点进行泰勒展开，得：

（紧接下屏）第六十三页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-64哈明(Hamming)公式（续1）令hk(k=0,1,…4)的系数为0,

得关于j、j的方程组:方程组有六个未知数，而只有五个方程，因此有无穷多解，取1为自由参量，则有:第六十四页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-65哈明(Hamming)公式（续2）当1=0时，得公式：

这是一个四阶隐式公式。完全类似，可导出如下形式的四阶显式公式：其中最常用的是：第六十五页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-66Hamming公式（预测一校正公式）

用显式公式（8-34）与隐式公式（8-33）联合可组成预测一校正公式

此式称为哈明（Hamming）公式，大量数值实验表明，此公式的数值稳定性在同类公式中几乎是最好的。为了进一步改善哈明公式计算结果的精度而又不致增加过多的计算量，可仿照带误差修正的阿当姆斯预测一校正公式，建立一个带误差修正的哈明公式：（紧接下屏）第六十六页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-67

前面介绍的四套预测一校正公式（8-28）、（8-29）、（8-32）、（8-33）都是实际计算中常用的方法，它们的共同特点是计算量省而精度高，但都必须用别的单步法为其提供出发值，由于解题过程是递推进行的，所以出发值的精度对以后的计算影响很大，故一般采用四阶龙格一库塔法求其出发值。Hamming公式（预测一校正公式）（续）

第六十七页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-68§4一阶方程组与高阶方程初值问题4.1一阶方程组

下面以两个方程的情形为例，介绍一阶方程组的数值解法。设初值问题：若采用向量记号，记：则初值问题（8-34）可表示为：

它与初值问题（8-1）形式上完全一致，前面介绍的解初值问题（8-1）的各种数值方法均可用于解一阶方程组初值问题（8-35），只要将公式中的y与f换成向量y与f即可。第六十八页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-69一阶方程组（续1）例如：解初值问题（8-35）的欧拉公式为：

即：写成分量形式即为：又如解初值问题（8-35）的四阶标准龙格-库塔公式为：第六十九页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-70一阶方程组（续2）写成分量形式为：

第七十页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-714.2高阶方程初值问题

高阶方程初值问题常转化为一阶方程组求解。例如：

引进新变量z=y，则此问题等价于

对其使用四阶标龙格-库塔法，计算公式为

第七十一页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-72高阶方程初值问题举例例5用四阶龙格-库塔方法在[0，1]上取步长h=0.2，求解二阶方程初值问题：解先将二阶方程化为一阶方程组。令z=y，则得方程组：使用四阶龙格-库塔公式，其相应的形式为：其中

第七十二页，共八十一页，2022年，8月28日第八章常微分方程数值解法8-73例5（续）表8-11

y-3y+2y=0解的近拟值xiyizik1l1k2l2k3l3k4l401.00001.00001.00001.00001.10001.10001.11001.11001.22201.22200.21.22141.22141.22141.22141.34351.34351.35581.35581.49261.49260.41.19181.19181.19181.19181.64101.64101.65591.65591.82301.82300.61