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

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

公式2.改进的欧拉公式3.龙格—库塔法4.亚当斯法5.算法的稳定性及收敛性2023/9/171§9.1引言包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y

及其各阶导数都是一次的,则称它是线性微分方程,否则称为非线性的。2023/9/172

在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解析解的方程是有限的，大多数的常微分方程不能给出解析解。譬如：

这个一阶微分方程就不能用初等函数来表达它的解。

2023/9/173再如，方程

虽然有解,仍需插值方法来计算各点的值2023/9/174从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题

(9.1)

在区间a≤x≤b

上的数值解法。

可以证明,如果函数在带形区域R:{a≤x≤b,-∞＜y＜∞}内连续，且关于y

满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使

对R内任意x及两个都成立,则方程(9.1)的解在

a,b

上存在且唯一。

2023/9/175数值方法的基本思想对常微分方程初值问题(9.1)式的数值解法，就是要算出精确解y(x)在区间

a,b

上的一系列离散节点处的函数值的近似值相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定h

为定数，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。

2023/9/176

对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，[1]它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息计算的递推公式。[2]建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等方法，对初值问题中的导数进行不同的离散化处理。

2023/9/177对于初值问题的数值解法，首先要解决的问题就是如何对微分方程中的导数进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算yi+1时只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格—库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法其代表是亚当斯法。

2023/9/178一阶微分方程的数值解可用数值微分、数值积分和泰勒展开等方法得到。

§9.2简单的数值方法与基本概念此公式称为欧拉公式2023/9/179

欧拉公式Euler

法的计算格式为：

2023/9/1710

Euler

格式用泰勒展开方法得到。2023/9/1711

Euler

格式用数值积分方法得到。

选择不同的计算方法计算上式的积分项就会得到不同的计算公式。

将方程的两端在区间上积分得，2023/9/1712

用左矩形方法计算积分项

代入(9.3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式

由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。

2023/9/1713Euler

公式的几何解释

欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)为通过点的一条积分曲线。积分曲线上每一点的切线的斜率等于函数在这点的值。

2023/9/1714Euler法的求解过程是:从初始点P0(即点(x0,y0))出发,作积分曲线y=y(x)

在P0点上切线(其斜率为

),与x=x1直线相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线方程为

当时,得这样就获得了P1点的坐标:

2023/9/1715同样,过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线的斜率直线方程为当时,得2023/9/1716当时,得由此获得了P2的坐标。重复以上过程,就可获得一系列的点:

。对已求得点以为斜率作直线

取2023/9/1717从图形上看,就获得了一条近似于曲线y=y(x)的折线。这样,从x0出发逐个算出2023/9/1718通常取(常数),则

微分方程的初值问题Euler

法的计算格式为：

用折线近似于曲线得到的计算公式称为欧拉公式2023/9/1719例9.1用欧拉法解初值问题

取步长

h=0.2,计算过程保留4位小数

解:欧拉迭代格式2023/9/1720例9.1用欧拉法解初值问题

取步长h=0.2,计算过程保留4位小数。

解:当k=0,x1=0.2时，已知x0=0，y0=1，有

y(0.2)

y1=0.2×1(4－0×1)＝0.8当k=1,x2=0.4时，已知x1=0.2，y1=0.8，有

y(0.4)

y2

=0.2×0.8×(4－0.2×0.8)＝0.6144当k=2,x3=0.6时，已知x2=0.4，y2=0.6144，有

y(0.6)

y3=0.2×0.6144×(4-0.4×0.6144)=0.4613

2023/9/1721clear;y=1,x=0,%初始化forn=1:10y=1.1*y-0.2*x/y,x=x+0.1,endy=1x=0y=1.1000x=0.1000y=1.1918x=0.2000y=1.2774x=0.3000y=1.3582x=0.4000y=1.4351x=0.5000y=1.5090x=0.6000y=1.5803x=0.7000y=1.6498x=0.8000y=1.7178x=0.9000y=1.7848x=1.0000解：计算结果如下：2023/9/1722对方程的两端在区间上积分得，代入(9.2)式,即可得到梯形公式为了提高精度,改用梯形方法计算其积分项，即

9.2.2梯形公式

由于欧拉（Euler）公式使用的是数值积分的左矩形方法精度很低，所以欧拉公式很粗糙。

2023/9/1723由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（9.5）是比欧拉公式(9.2)精度高的一个数值方法。9.2.2梯形公式梯形公式2023/9/1724(9.5)

(9.2)和(9.5)式粗看差不多，

(9.2)

Euler

公式梯形公式但(9.5)式的右端也含有未知的yi+1,它是一个关于yi+1的方程,这类数值方法称为隐式方法。相反地,欧拉公式是关于yi+1

的一个直接的计算公式，这类数值方法称为显式方法。隐式公式,需用迭代法求解。2023/9/17259.2.3两步欧拉公式

对方程的两端在区间上积分得

(9.6)

改用中矩形公式计算其积分项，即

代入上式,并用yi

近似代替式中y(xi)即可得到两步的欧拉公式2023/9/1726

前面介绍的欧拉方法、梯形方法，它们都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;

无论是单步法,还是两步法只要是直接法就需估计它的误差。

而中矩形公式(9.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息,故称其为两步欧拉公式。

2023/9/1727定义9.1在yi准确的前提下,即时,用数值方法计算yi+1的误差,称为该数值方法计算时yi+1

的局部截断误差。9.2.4．欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度,因此引入局部截断误差和收敛阶的概念。在假设前i-1步精确的前提下，得到的第i步的误差称为局部截断误差2023/9/1728假定，由欧拉公式则有因此有

局部截断误差2023/9/1729定义9.2如果数值方法的局部截断误差为,则称这种数值方法的阶数是P。显然，步长(h<1)越小,P越高,则局部截断误差越小,计算精度越高。两步欧拉公式(中矩形)欧拉公式的局部截断误差为,欧拉方法仅为一阶方法。是比欧拉公式精度高的一个数值方法。2023/9/1730设

前两步准确,则两步欧拉公式

把y(xi-1)在xi处展开成Taylor级数，即

①把y(xi-1)代人①，即

(9.8)2023/9/1731再将精确值y(xi+1)在xi处进行泰勒展开(9.8)(9.9)所以,由(9.8)和(9.9)可得两步欧拉公式的局部截断误差为：

两步欧拉公式为二阶方法，故比一步欧拉公式精度高。2023/9/1732显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量非常大。将欧拉公式和梯形公式综合，便可得到改进的欧拉公式。9.2.5改进的欧拉公式欧拉公式梯形公式2023/9/1733先用欧拉公式(9.2)求出一个初步的近似值,称为预测值,它的精度不高,再用梯形公式对它校正一次,即再迭代一次,求得yi+1,称为校正值,欧拉公式梯形公式这种预测-校正方法称为改进的欧拉公式。2023/9/1734(9.10)

预测

校正左矩形公式梯形公式改进的欧拉公式右矩形公式一次左矩形公式，一次梯形公式2023/9/1735可以证明,改进的欧拉公式(9.10)的精度为二阶。这是一种一步显式格式,它可以表示为嵌套形式。或者表示成下列平均化形式2023/9/1736

改进的欧拉公式左矩形公式右矩形公式左、右矩形公式的平均值，即梯形公式，又称为改进的欧拉公式。2023/9/17379.2.6改进欧拉法算法实现（1）计算步骤①输入,h,N②使用以下改进欧拉法公式进行计算

③

输出，并使转到

②

直至n>N结束。

2023/9/1738（2）改进欧拉法的流程图

2023/9/1739例9.2用改进欧拉法解初值问题区间为

0,1,取步长h=0.1

解:改进欧拉法的具体形式本题的精确解为2023/9/1740clearx=0,yn=1%初始化forn=1:10yp=yn+0.1*(yn-2*x/yn);%预测x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end2023/9/1741例9.3对初值问题

求得的近似解为

并证明当步长h

0时,yn收敛于精确解整理成显式

反复迭代,得到证明:解初值问题的梯形公式为证明用梯形公式2023/9/1742反复迭代,得到由于，有

证毕证毕2023/9/1743§9.3龙格-库塔（Runge-Kutta）法只有对平均高度f(x,y)提供一种算法，便可得到一种计算格式。9.3.1龙格-库塔法的基本思想2023/9/1744Euler

公式可改写成改进的Euler公式又可改写成用左端点的高度作为平均高度得到的计算格式用左、右端点高度的平均值作为平均高度得到的计算格式2023/9/1745上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值(高度)的线性组合得出y(xi+1)的近似值yi+1,而且增加计算f(x,y)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次f(x,y)的值,为一阶方法局部截断误差为。改进欧拉公式需计算两次f(x,y)

的值，它是二阶方法。它的局部截断误差为。2023/9/1746

于是可考虑在内多预报几个点的函数值(高度)，然后将其加权平均作为平均高度，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi

处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。则可构造出更高精度的计算格式，这就是龙格—库塔（Runge-Kutta）法的基本思想。2023/9/1747只有多计算几个点的函数值，平均高度就更精确.龙格-库塔法的基本思想或多预报几个点的斜率值，将其加权平均作为平均斜率，平均斜率就更精确。2023/9/17489.3.2二阶龙格—库塔法式中:比照改进的欧拉法2023/9/1749式中K可看作是y=y(x)在区间上的平均斜率。所以可得计算公式为：（9.14）

也即（9.13）其中：由微分中值定理，存在点，使得2023/9/1750二元函数的Taylor展开式2023/9/1751将以上结果代入将在x=xi

处进行二阶Taylor展开：

（9.15）

2023/9/1752对式(9.15)和(9.16)进行比较系数后可知,只要成立,格式(9.14)的局部截断误差就等于式(9.17)中具有三个未知量,但只有两个方程,因而有无穷多解。2023/9/1753式(9.17)中，若取,则P=1，这是无穷多解中的一个解，将以上所解的值代入式(9.14)并改写可得

不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（9.17）有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。2023/9/1754若取,则，此时二阶龙格-库塔法的计算公式为

此计算公式称为变形的二阶龙格—库塔法。式中为区间的中点。2023/9/17559.3.3三阶龙格-库塔法

为了进一步提高精度，设除外再增加一点并用三个点的斜率k1,k2,k3

加权平均得出平均斜率k*的近似值，这时计算格式具有形式:（9.18）

为了预报点的斜率值k3,在区间内有两个斜率值k1和k2

可以用,可将k1,k2加权平均得出上的平均斜率,从而得到的预报值2023/9/1756于是可得

运用Taylor展开方法选择参数,可以使格式(9.18)的局部截断误差为,即具有三阶精度.2023/9/1757局部截断误差为,即具有三阶精度，这类格式统称为三阶龙格—库塔方法。（9.19）

是其中的一种，称为库塔（Kutta）公式。

2023/9/17589.3.4四阶龙格—库塔法

如果需要再提高精度，用类似上述的处理方法，只需在区间上用四个点处的斜率加权平均作为平均斜率k*的近似值，构成一系列四阶龙格—库塔公式。具有四阶精度，即局部截断误差是。由于推导复杂，这里从略，只介绍最常用的一种四阶经典R—K公式。

2023/9/1759

四阶经典R—K公式

（9.20）

2023/9/17609.3.5四阶龙格—库塔法算法实现(1)

计算步骤①输入,h,N②使用龙格—库塔公式（9.20）计算出y1③输出，并使转到②直至n>N结束。2023/9/1761(2)四阶龙格—库塔算法流程图2023/9/1762例9.4取步长h=0.2，用经典R-C格式求解初值问题解:由四阶龙格-库塔公式可得2023/9/1763可同样进行其余yi

的计算。本例方程的解为，从表中看到所求的数值解具有4位有效数字。

龙格—库塔法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用龙格—库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。

2023/9/1764

以经典的四阶龙格-库塔法(9.20)为例。从节点xi出发，先以h为步长求出一个近似值，记为，由于局部截断误差为，故有

当h值不大时，式中的系数c可近似地看作为常数。9.3.6变步长的龙格-库塔法在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。2023/9/1765然后将步长折半,即以为步长,从节点xi出发,跨两步到节点xi+1,再求得一个近似值,每跨一步的截断误差是,因此有这样由此可得这表明以作为的近似值，其误差可用步长折半前后两次计算结果的偏差来判断所选步长是否适当2023/9/1766当要求的数值精度为ε时：（1）如果Δ>ε，反复将步长折半进行计算，直至Δ<ε为止,并取其最后一次步长的计算结果作为（2）如果Δ<ε，反复将步长加倍，直到Δ>ε为止，并以上一次步长的计算结果作为。这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要反复判断Δ，增加了计算工作量，但在方程的解y(x)变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。2023/9/1767§9.4算法的稳定性及收敛性9.4.1稳定性稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。

2023/9/1768

当在某节点上xi的yi值有大小为δ的扰动时，如果在其后的各节点上的值产生的偏差都不大于δ，则称这种方法是稳定的。

稳定性不仅与算法有关，而且与方程中函数f(x,y)也有关，讨论起来比较复杂。为简单起见，通常只针对线性模型方程来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求解。2023/9/1769先考察显式Euler方法的稳定性。模型方程的Euler公式为将上式反复递推后，可得或式中2023/9/1770

要使yi

有界，其充要条件为即由于，故有（9.26）可见，如欲保证算法的稳定，显式Euler格式的步长h的选取要受到式（9.26）的限制。的绝对值越大，则限制的h值就越小。用隐式Euler格式右矩形公式，对模型方程的计算公式为可化为2023/9/1771由于,则恒有，故恒有因此，隐式Euler

格式是绝对稳定的（无条件稳定）的（对任何h

>

0）。2023/9/1772常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值yi来近似替代y(xi)，这种近似替代是否合理，还须看分割区间的长度h越来越小时，即时，是否成立。若成立，则称该方法是收敛的，否则称为不收敛。9.4.2收敛性这里仍以Euler方法为例，来分析其收敛性。Euler格式2023/9/1773设表示取时,按Euler公式的计算结果,即Euler方法局部截断误差为：设有常数，则（9.27）

总体截断误差（9.28）

又由于f(x,y)关于y满足李普希茨条件,即2023/9/1774代入(9.28)上式，有再利用式（9.27），式（9.28）即上式反复递推后，可得（9.29）

2023/9/1775设（T为常数）

因为

所以把上式代入式（9.29），得若不计初值误差，即，则有（9.30）

式(9.30)说明,当时,,即,所以Euler方法是收敛的，且其收敛速度为，即具有一阶收敛速度。同时还说明Euler方法的整体截断误差为，因此算法的精度为一阶。2023/9/1776§9.5亚当姆斯

Adams

方法龙格-库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率值，计算量比较大。考虑到计算yi+1

之前已得出一系列节点上的斜率值，能否利用这些已知的斜率值来减少计算量呢？这就是亚当姆斯（Adams）方法的设计思想。

9.5.1亚当姆斯格式2023/9/1777（9.21）

选取参数λ，使格式（9.21）具有二阶精度。精确值近似值前三项相同2023/9/1778相比较,需取格式(9.21)具有二阶精度。称之为二阶亚当姆斯格式。这样导出的计算格式2023/9/1779（9.22）

类似地可以导出三阶亚当姆斯格式。

2023/9/1780类似地四阶亚当姆斯格式。

(9.22)上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点的斜率值来预报区间上的平均斜率是个外推过程，效果不够理想。为了进一步改善精度，变外推为内插，即增加节点xi+1的斜率值来得出上的平均斜率。譬如考察形如2023/9/1781（9.23）

(9.23)为隐式格式,可见要使格式(9.23)具有二阶精度,需令,这样就可构造二阶隐式亚当姆斯格式①代人①得2023/9/1782（9.23）

二阶隐式亚当姆斯格式类似可导出三阶隐式亚当姆斯格式其实是梯形格式。和四阶隐式亚当姆斯格式

（9.24）

2023/9/17839.5.2亚当姆斯预报-校正格式

参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式(9.22)和隐式(9.24)相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报-校正格式（9.25）

预报：

校正：

2023/9/1784这种预报-校正格式是四步法，它在计算yi+1时不但用到前一步的信息，而且要用到再前面三步的信息，因此它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格—库塔法提供开始值。

2023/9/1785例9.5取步长h=0.1,用亚当姆斯预报-校正公式求解初值问题

的数值解。解:

用四阶龙格-库塔公式求出初值,计算得：表中的,yi和y(xi)分别为预报值、校正值和准确解(),以比较计算结果的精度。再使用亚当姆斯预报-校正公式(9.25),2023/9/1786§9.6一阶常微分方程组与高阶方程我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。

9.6.1一阶常微分方程组对于一阶常微分方程组的初值问题

（9.31）

可以把单个方程中的f和y看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。2023/9/1787设为节点上的近似解，则有改进的Euler格式为预报：校正：

（9.32）

又，相应的四阶龙格—库塔格式（经典格式）为（9.33）

2023/9/1788式中

（9.34）

2023/9/1789把节点xi上的yi和zi值代入式(7.34),依次算出

,然后把它们代入式(7.33),算出节点xi+1上的yi+1

和zi+1值。对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初