《计算方法》课件第7章

第7章常微分方程的数值解法7.1引言7.2欧拉方法及改进的欧拉方法7.3龙格－库塔方法7.4线性多步法 7.1引言

函数是客观事物的内部联系在数量方面的反映，利用函数关系式可以对客观事物的规律性进行研究，因此寻求函数关系式在实践中具有重要意义。然而工程实践和科学研究中，许多问题仅能够提供包含有要找的函数及其导数的关系式，也就是微分方程。通过微分方程的求解可以得到函数关系式，这在高等数学中已经给出了一些典型微分方程求解的基本方法。但实际中遇到的微分方程往往比较复杂，在很多情况下，不能给出解的解析表达式。有时候即使能用解析表达式表示，又因为计算量太大而不实用。更重要的是，许多实际问题，并不需要方程解的表达式，而仅仅需要得到一些点上的近似值即可。因此，有必要学习和研究微分方程的数值解法。在微分方程中，如果未知函数是一元函数，则微分方程称做常微分方程，未知函数是多元函数的称为偏微分方程。而微分方程中，出现的未知函数最高阶导数的阶数称为微分方程的阶。如果给定了微分方程，并给出了初始条件，则被称为微分方程初值问题;如果给出的是边界条件，则称之为边值问题。本章研究微分方程的数值解法，主要讨论微分方程中最简单的一类问题，即一阶常微分方程初值问题。一阶常微分方程初值问题描述如下:(7.1)假设式(7.1)在区间［a，b］上存在唯一且足够光滑的解y(x)。式(7.1)的数值求解就是寻求解y(x)在一系列离散点a≤x0<x1<…<xn≤b(也称为节点)处的近似函数值yi=f(xi)(i=0，1，2，…，n)。相邻两个节点间的距离——节点间距，也称为步长，用hi表示，即hi=xi+1－xi(i=0，1，2，…，n－1)。通常采用等距节点，即hi=h(常数)，此时可以记xi=a+ih(i=0，1，2，…，n)。7.2欧拉方法及改进的欧拉方法

7.2.1欧拉方法

1.欧拉公式

欧拉方法是解一阶常微分方程初值问题的基本方法，可以利用数值微分的概念来推导欧拉公式。

关于一阶微分方程初值问题，其目标是计算函数在x1，x2，…，xn点上的函数值y1，y2，…，yn。为此，利用数值微分中的向前差商近似导数的方法可得由此可得根据式(7.1)可以得到y′(x0)=f(x0，y0)及y(x0)=y(a)=y0，代入式(7.2)得并将其记为(7.2)(7.3)式中,x0、y0、h及f(x，y)的表达式都是已知的，所以可以根据式(7.3)计算出节点x1处的函数近似值y1。用同样的方法，可以得到节点x2处的函数近似值y2=y1+hf(x1，y1)，…，直到节点xn处的函数近似值yn=yn－1+hf(xn－1，yn－1)，则计算过程的通用表达式为式(7.4)即为求解一阶常微分方程初值问题的欧拉公式。

例7－1取步长h=0.02，用欧拉公式求一阶微分方程初值问题：解根据欧拉公式及已知条件可得由于h=0.02，区间为［0，0.1］，所以节点为x0=0，x1=0.02，x2=0.04，x3=0.06，x4=0.08，x5=0.1，计算过程略，计算结果如表7.1所示。表7.1例7－1的计算结果

2.欧拉公式的几何意义

设一阶常微分方程的解函数y(x)为如图7.1中所示的曲线，根据一阶常微分方程初值问题的已知条件，曲线y(x)上的初始点(x0，y0)已知，故过该初始点以f(x0，y0)为斜率,可以得到一条直线，且该直线的方程为

y=y0+f(x0，y0)(x－x0)

该直线与直线x=x1的交点设为P1，则P1的纵坐标为记为y1，即此与利用欧拉公式得到的值相同。然后过点(x1，y1)以f(x1，y1)为斜率，可以得到一条直线，且该直线方程为同样可以得到该直线与直线x=x2的交点P2的纵坐标，记为y2，即此与利用欧拉公式得到的值也相同。依此类推，得到图7.1欧拉公式的几何意义

3.欧拉方法的局部截断误差

定义7.1假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri+1=y(xi+1)－yi+1，称为局部截断误差，可用图7.2表示，图中|APi+1|即为欧拉方法在xi+1点的局部截断误差。图7.2欧拉方法的局部截断误差示意图

定义7.2如果数值方法得到的局部截断误差为O(hp+1)，则称该方法具有p阶精度或称为p阶的。

根据局部截断误差的定义，可以计算欧拉公式的局部截断误差。将y(xi+1)在xi处泰勒展开可得(7.5)由欧拉公式可得(7.6)根据局部截断误差的定义有则由式(7.5)和式(7.6)可得(7.7)式(7.7)即为欧拉方法的局部截断误差，故其具有一阶精度。7.2.2改进的欧拉方法

从欧拉方法的局部截断误差和欧拉方法的几何意义都可以看出，欧拉方法的误差较大，为此研究其改进方法，以获得精度较好的计算方法。

一阶微分方程初值问题的微分方程为对其在小区间［xi，xi+1］(i=0，1，…，n－1)上进行积分，可得即(7.8)根据问题的初始条件，如果能够计算出式(7.8)右端的积分值，则可以获得微分方程的数值解,故求解一阶微分方程初值问题变成了求解函数积分的问题。

1.应用矩形公式

对式(7.8)右端的积分应用左矩形公式可得将其代入式(7.8)可得将y(xi)记做yi，y(xi+1)记做yi+1，则有可以看出式(7.9)与欧拉公式完全相同，这说明用数值积分的方法也可以推导得到欧拉公式。

2.应用梯形公式

对式(7.8)右端的积分应用梯形公式可得将其代入式(7.8)可得将y(xi)记做yi，y(xi+1)记做yi+1，则有(7.10)式(7.10)是应用了数值积分中的梯形公式推导得到的，因此也被称为梯形公式。根据局部截断误差的定义可以得到梯形公式的局部截断误差为(7.11)同样，将y(xi+1)在xi处泰勒展开可得(7.12)其中，将式(7.12)代入式(7.11)有(7.13)将y′(xi+1)在xi处泰勒展开可得(7.14)其中,xi<τi<xi+1。将式(7.14)代入式(7.13)得(7.15)

3.改进的欧拉公式

欧拉公式计算简单，但误差较大。梯形公式虽可以提高精度，但为隐式公式，缺点也比较明显(计算量增大了)。为此，将欧拉公式和梯形公式联合使用，克服梯形公式的缺点，由此得到了改进的欧拉公式。改进的欧拉公式就是先用欧拉公式得出一个y(xi+1)的近似值，用yi+1来表示，称为预估值。然后对预估值yi+1使用梯形公式进行调整，得到较为精确的近似值yi+1，称之为校正值。整个计算公式为(7.16)式(7.16)中第一式为预估算式，第二式为校正算式，整个公式称为改进的欧拉公式。为了便于编写程序，常将式(7.16)改写为(7.17)

例7－2取步长h=0.1，用欧拉方法和改进的欧拉方法求如下一阶微分方程初值问题：

解由于h=0.1，区间为［0，1.5］，x0=0，所以需计算函数值的节点为xk=0+0.1k(k=1，2，…，15)。

(1)用欧拉方法求解，其计算公式为

(2)用改进的欧拉方法求解，其计算公式为将已知数据代入上面的计算公式，可以得到微分方程的数值解，计算过程略，计算结果见表7.2。 由于该微分方程的解析解为 ，为了比较算法计算结果的误差，将真实值y(xi)(i=0，1，2，…，15)也列入表7.2中。从表中数据可以看出，改进的欧拉方法的计算结果更接近真实值。表7.2例7－2的计算结果

7.3龙格－库塔方法

7.3.1龙格－库塔方法的基本思想

对微分方程y′=f(x，y)，在小区间［xi，xi+1］上应用拉格朗日中值定理，可得(7.18)令ξ=xi+θh，θ∈(0，1)，则式(7.18)变为(7.19)根据式(7.19)，如果能够计算出f(xi+θh，y(xi+θh))，则可以根据一阶微分方程初值问题的条件递推得到问题的数值解，而f(xi+θh，y(xi+θh))在几何意义上表示曲线上点的切线斜率，因此令K*=f(xi+θh，y(xi+θh))(称其为区间［xi，xi+1］上的平均斜率)，则式(7.19)变为(7.20)因此只要对K*提供一种算法，就可以求得微分方程的数值解。下面根据该观点对欧拉方法及改进的欧拉方法进行分析。(7.21)将y(xi)记做yi，y(xi+1)记做yi+1，则有(7.22)该式即欧拉公式，或者说欧拉公式以点(xi，yi)上的斜率来计算K*。当取θ=0和θ=1时，可以得到两个点(xi，yi)和(xi+1，y

i+1)上的斜率，并分别令其为K1和K2，则有(7.23)其中,点(xi+1，yi+1)处的斜率K2是由已知信息预测得到的。 如果取点(xi，yi)和点(xi+1，yi+1)上的斜率的算术平均值作为区间［xi，xi+1］上的平均斜率，即(7.24)则由式(7.20)、式(7.23)和式(7.24)，并将y(xi)和y(xi+1)分别记作yi和yi+1，得(7.25)

这显然是改进的欧拉公式，或者说改进的欧拉公式以两个点(xi，yi)和(xi+1，yi+1)上的斜率的算术平均值来计算K*。

已经知道改进的欧拉公式具有二阶精度，而欧拉公式只有一阶精度，结合上面的分析可知，用两个点斜率的算术平均作为平均斜率比只取一个点的斜率作为平均斜率的计算精度要高。由此，如果设法在区间［xi，xi+1］内多预测几个点的斜率值，然后将它们加权平均作为平均斜率K*，则有可能构造出具有更高精度的计算公式，这就是龙格－库塔方法的基本思想。7.3.2二阶龙格－库塔方法

在区间［xi，xi+1］上，改进的欧拉方法使用了xi、xi+1两个节点上的斜率值。现将其推广到一般情况，取区间［xi，xi+1］上的点xi、xi+p，其中xi+p=xi+ph(0<p≤1)。设xi

和xi+p两个点上的斜率分别为K1和K2，取其加权平均作为平均斜率K*的近似值，即(7.26)其中λ1和λ2为权值。显然K1=f(xi，yi)，而K2=f(xi+p，yi+p)，对于其中的yi+p，应用欧拉公式计算得(7.27)由式(7.20)、式(7.26)和式(7.27)可得(7.28)式(7.28)中，λ1、λ2、p为待定的未知数，如果能够确定λ1、λ2、p的值，就可以计算微分方程的数值解。为了确定未知数的值，可以利用式(7.28)具有二阶精度的条件，为此，首先计算式(7.28)的局部截断误差。根据泰勒公式有根据式(7.28)有(7.30)对进行泰勒展开可得将其代入式(7.30)得(7.31)由式(7.29)和式(7.31)可得式(7.28)的局部截断误差为(7.32)式(7.28)具有二阶精度，则其局部截断误差Ri+1=O(h3)，所以根据式(7.32)有(7.33)方程组(7.33)为三个未知数，两个方程，因此方程组有无穷多个解，每一个解都可以得到一个具有二阶精度的龙格－库塔公式。特别地，当λ1=λ2=1/2，p=1时,式(7.28)可写为(7.34)式(7.34)显然是改进的欧拉公式。此外，当p=1/2时，可以得到λ1=0，λ2=1/2，将其代入式(7.28)有(7.35)式(7.35)称为变形的欧拉公式。7.3.3高阶龙格－库塔方法

1.三阶龙格－库塔方法

根据龙格－库塔方法的思想，为了进一步提高微分方程数值求解的精度，可以在区间［xi，xi+1］上，除了点xi、xi+p外，再增加节点。如再增加一个点xi+m，xi+m=xi+mh，其中p<m≤1，则利用该三点处的斜率K1、K2、K3的加权平均值作为K*近似值，即(7.36)可得计算公式为(7.37)式中K1和K2的取法与二阶龙格－库塔方法相同，分别为而K3=f(xi+m，yi+m)需要得到点xi+m处的函数值y(xi+m)，为此用二阶龙格－库塔公式对yi+m进行估计,即(7.38)因此可得(7.39)由此得计算公式为(7.40)式(7.40)中，λ1、λ2、λ3、p、m、μ1、μ2为待定的未知数，如果能够确定它们的值，就可以得到微分方程的数值解。为了确定未知数的值，与二阶龙格－库塔公式的推导过程类似，利用式(7.40)具有三阶精度的条件，借助泰勒展开，计算其局部截断误差，从而得到待定未知数满足的条件(方程)。具体推导过程略，得到的待定未知数满足的条件为(7.41)显然,由式(7.41)构成的方程组有无穷多个解，每一个解都可以得到一个具有三阶精度的龙格－库塔公式，因此称满足条件(7.41)的计算公式(7.40)为三阶龙格－库塔公式。常用的三阶龙格－库塔公式为(7.42)

2.四阶龙格－库塔方法

与推导三阶龙格－库塔公式的方法类似,可以得到四阶龙格－库塔公式，只需在区间［xi，xi+1］上用四个点处的斜率的加权平均作为K*的近似值即可，计算公式的推导过程略。常用的经典四阶龙格－库塔公式如下：(7.43)基于龙格－库塔公式的思想，也可推导出其更高阶的公式，龙格－库塔公式的一般表达式为(7.44)其中，需要说明的是，龙格－库塔公式的推导基于泰勒展开方法，所以它要求解函数有很好的光滑性，即y(x)要具有所要求的导数。不然，高阶龙格－库塔公式可能不如低阶龙格－库塔公式的效果好。通常来说，对于一般实际问题，四阶龙格－库塔方法即可达到一般的精度要求。

例7－3分别用三阶和四阶龙格－库塔公式计算例7－2的初值问题。

解由于h=0.1，区间为［0，1.5］，x0=0，所以需计算函数值的节点为xi=0+0.1i(i=1，2，…，15)。

(1)采用三阶龙格－库塔公式计算。

由公式(7.42)可得(7.43)基于龙格－库塔公式的思想，也可推导出其更高阶的公式，龙格－库塔公式的一般表达式为(7.44)其中，需要说明的是，龙格－库塔公式的推导基于泰勒展开方法，所以它要求解函数有很好的光滑性，即y(x)要具有所要求的导数。不然，高阶龙格－库塔公式可能不如低阶龙格－库塔公式的效果好。通常来说，对于一般实际问题，四阶龙格－库塔方法即可达到一般的精度要求。

例7－3分别用三阶和四阶龙格－库塔公式计算例7－2的初值问题。

解由于h=0.1，区间为［0，1.5］，x0=0，所以需计算函数值的节点为xi=0+0.1i(i=1，2，…，15)。

(1)采用三阶龙格－库塔公式计算。

由公式(7.42)可得

(2)采用四阶龙格－库塔公式计算。

由公式(7.43)可得其中,i=0，1，2，…，14。将x0=0，y0=1代入进行计算，计算过程略，计算结果见表7.3。同样，根据该微分方程的解析解，将真实值y(xi)(i=0，1，2，…，15)也列于表7.3中。

从表7.3可以看出，四阶龙格－库塔公式的计算结果更接近于真实值。表7.3例7－3的计算结果 7.4线性多步法

7.4.1线性多步法的基本思想

前面所讨论的求解一阶微分方程初值问题的欧拉方法及龙格－库塔方法，在计算节点xi+1上的函数值yi+1时，仅利用了节点xi上的函数值yi，而此时已经计算得到了节点xi+1前面的一系列点上的函数值的近似值y0，y1，…，yi，但这些近似值没有被利用，这样的方法被称为单步法。相反，在计算节点xi+1上的函数值yi+1时，不仅利用了节点xi上的函数值yi，还利用了节点xi前面的一系列节点上的函数值的近似值yi－1，yi－2，…，这样的计算方法称为多步法。例如r步法，即计算节点xi+1上的函数值yi+1时，用到了yi，yi－1，yi－2，…,yi－r+1。

多步法计算中，由于使用了前面多步的信息来计算当前节点的函数值，可以获得比单步法更高的精度。最常用的多步法是线性多步法。r步线性多步法的公式可表示为(7.45)(7.46)在已知节点x0，x1，…，xn及y0，y1，…，yi的情况下，选取不同的插值节点做插值多项式就会得到不同的数值解法，常用的阿达姆斯方法就是基于该思想构造的。7.4.2阿达姆斯显式公式

针对式(7.46)在已计算出y0，y1，…，yi的情况下，选取r个节点xi，xi－1，…，xi－r+1进行拉格朗日插值，可得(7.47)其中，将式(7.47)作为插值函数P(x)代入式(7.46)，并记yk=y(xk)，可得(7.48)令则式（7.48）可表示为：(7.49)式(7.49)中，当r=4时可以得到更加具体的四步法公式：(7.50)其中,由于是等间距节点，设x=xi+th(x∈［xi，xi+1］，t∈［0，1］)，将其代入上式可得同理可得则式(7.50)的具体表达式为(7.51)式(7.51)即为常用的四步阿达姆斯显式公式，它具有四阶精度，其局部截断误差可以根据拉格朗日插值余项计算得到,即(7.52)7.4.3阿达姆斯隐式公式

针对式(7.46)在已计算出y0，y1，…，yi的情况下，选取r个节点xi+1，xi，…，xi－r+2进行拉格朗日插值，可得(7.53)其中，将式(7.53)作为插值函数P(x)代入式(7.46)，并记yk=f(xk)，可得(7.54)令，则式(7.54)可表示为(7.55)很显然,式