《计算方法》课件1第8章

第八章常微分方程初值问题的数值解法8.1引言8.2欧拉方法及其改进8.3龙格—库塔方法8.4线性多步法8.5算法的稳定性及收敛性8.6一阶常微分方程组与高阶方程8.7解微分方程的波形松弛方法8.8微分方程边值问题的数值方法习题88.1引言

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

解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解的。另外，在许多实际问题中，并不需要方程解的表达式，而仅仅需要获得解在若干点上的近似值即可。因此，研究常微分方程的数值解法很有必要。在本章里，将向大家着重介绍一阶常微分方程初值问题：

(8.1)

在区间［a,b］上的数值解法。

下面介绍一些与数值解法有关的内容。

定理8.1如果函数f(x，y)在带形区域D={(x,y)|a≤x≤b，

y∈R}上有定义且连续，并且满足李普希兹条件，即当

(x,

y1)，(x，y2)∈D时，有

|f(x，y1)－f(x，y2)|≤L|y1－y2|

其中，L(0<L<+∞)称为李普希兹常数(它与x,y无关)，则对任何(x0,y0)∈D，初值问题((8.1)式)在［a,b］上存在唯一连续可微解y=y(x)。

定理8.2如果函数f(x，y)在带形区域D={(x，y)|a

≤x≤b，y∈R}上满足李普希兹条件，则初值问题((8.1)式)的数值解是稳定的。

(证略。)

所谓常微分方程初值问题的数值解法，就是要算出精确解y=y(x)在区间［a,b］上的一系列离散节点a≤x0<x1<x2<…

<xn<…≤b处的函数值y(x0),y(x1),…，y(xn)的近似值y0，y1，…，yn，…。相邻两个节点的间距Δxi=xi－xi－1称为步长，步长可以相等，也可以不等。在本章中，一般取等步长并用字母h表示，并且记f(xn，y(xn))=y′(xn)。求初值问题的数值解法采用的是“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，算出yn后再算yn+1。这种数值方法有单步法和多步法之分。单步法是在计算yn+1

时，仅利用yn；多步法是在计算yn+1时不仅要用到yn，还要用到yn－1，yn－2,…，一般k步法要用到yn,yn－1,…，

yn－k+1

。单步法的代表方法是龙格—库塔方法；多步法的代表是亚当斯方法。

单步法和多步法又有显式和隐式之分。显式单步法形如yn+1=yn+hφ(xn,yn,h)；隐式单步法形如yn+1=yn+hφ(xn,yn,

yn+1,h)。对多步法来说，显式与隐式方法与单步法相似。

8.2欧拉方法及其改进

8.2.1欧拉公式

欧拉(Euler)方法是解初值问题((8.1)式)最简单的数值方法。推导欧拉公式有多种方法，譬如利用差商、数值积分、数值微分、泰勒展开法等。下面我们以数值积分为例进行推导。对(8.1)式的第一式y′=f(x,y)两端在区间［xn,xn+1］上进行积分，得

即

(8.2)选择不同的计算方法计算积分项,

就会得到一系列不同形式的欧拉公式。

1.向前欧拉法

用左矩形公式计算积分项，得

代入(8.2)式中，并用yn近似代替式中y(xn),即可得到向前欧拉公式:

(8.3)

这是一种显式形式的方程,因此也称为显式欧拉公式。

2.向后欧拉法

用右矩形公式计算积分项，得

代入(8.2)式中，并用yn+1近似代替式中y(xn+1)，即可得向后欧拉公式:

yn+1=yn+hf(xn+1,yn+1)

(8.4)

这是一种隐式形式的方程,因此也称为隐式欧拉公式。由于数值积分的矩形方法精度很低，因此欧拉公式比较粗糙。

欧拉方法的几何意义十分清楚，如图8-1所示。(8.1)式的解曲线y=y(x)过P0(x0,y0)点，从P0出发以f(x0,y0)为斜率做直线段，与x=x1相交于P1(x1，y1)，显然有y1=y0+hf(x0,y0)。同理，再从P1出发，以f(x1,y1)为斜率做直线段，与x=x2相交于P2(x2,y2)。依次类推，可得一条折线P0P1…Pn…，作为解曲线y=y(x)的近似曲线，故欧拉方法又称为欧拉折线法。图8-1

3.梯形方法

为了提高精度，用梯形公式计算积分项,

得

代入(8.2)式中，并用yn近似代替式中y(xn)，yn+1近似代替式中y(xn+1)，即可得到梯形公式:

(8.5)由于数值积分的梯形公式要比矩形公式精度高，因此梯形公式((8.5)式)要比欧拉公式((8.3)式)精度高。梯形法是一种隐式单步法，从n=0开始，每步都要解关于yn+1的一个方程。一般来说，这是一个非线性方程，因此要用迭代法求解。

4.改进欧拉方法

显式欧拉公式计算工作量小，但精度低。梯形公式虽提高了精度，但为隐式公式，需用迭代法求解，计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。

具体地说，就是先用欧拉公式求出一个初步的近似值

y［0］n+1，称之为预估值。预估值的精度可能很差，再用梯形公式将它校正一次，即迭代一次，求得的yn+1称为校正值，由这种预测—校正方法得到的公式称为改进的欧拉公式：

(8.6)其中，第一式称为预估算式；第二式称为校正算式。有时为了计算方便，将(8.6)式改写为按照上述思想，如果预估时用其他公式(譬如中矩形公式、辛普森公式等)，则可以得到其他形式的预估—校正式。例如，若采用中矩形公式，则有下述公式：8.2.2单步法的局部截断误差和阶

定义8.1称

en=y(xn)－yn

为某一数值方法在点xn处的整体截断误差；称

Rn,h=y(xn+1)－y(xn)－hφ(xn,y(xn),h)

为显式单步法式yn+1=yn+hφ(xn,yn,h)在xn+1处的局部截断误差。其中，y(x)为初值问题的精确解。

Rn，h之所以称为局部的，是因为如果假设y(xn)=yn，即第n步以及以前各步都没有误差，则由式yn+1=yn+hφ(xn,yn,h)所得的yn+1与y(xn+1)之差为

即在假定的y(xn)=yn条件下，Rn，h=y(xn+1)－yn+1，这就是

Rn,h称为局部的含义。

定义8.2若数值方法的局部截断误差为O(hp+1)，则称这种方法为p阶的。

通常p越大，h越小，则截断误差越小，数值方法越精确。

可以证明，欧拉方法是一阶方法，改进的欧拉方法以及梯形方法是二阶方法。下面仅给出欧拉方法的证明。

设y(xn)=yn，把y(xn+1)在xn处展开成泰勒公式，即则由欧拉公式可得

两式相减得欧拉公式的局部截断误差为若y(x)在［a,b］上充分光滑，且令,

则

故欧拉方法是一阶方法。

例8.1应用向前欧拉公式求初值问题:

取步长h=0.1，将计算结果与精确解y=x+e-x对照。

解将区间［0,1］进行10等分，h=0.1，xn=nh(n=0.1,

…，10)。

向前欧拉公式为

数值解yn与精确解y(xn)及误差列于表8.1。

例8.2用改进欧拉法求初值问题：

取步长h=0.1，将计算结果与精确解y=e-x比较。(0≤x≤1)

解将区间［0,1］进行10等分，步长h=0.1，节点xn=nh(n=0.1,…,10)。

应用改进的欧拉公式，得

由此，y1=0.905,y2=0.9052,…，y10=0.90510。

数值解yn与精确解y(xn)及误差列于表8.2。

例8.3有初值问题：

试证明：

(1)用梯形公式求得的数值解为

(2)当步长h→0时，yn收敛于精确解y=e-x。(0≤x≤1)

证明

(1)应用梯形公式:

整理上式，得由此公式递推可得

于是

(2)设区间是等距划分的，对于任意给定的节点x=xn=nh，步长。显然当h→0的同时，n→∞，

由此

证毕。8.3龙格—库塔方法

本节中将向大家介绍求解初值问题的一类高精度的单步法——龙格—库塔(Runge-Kutta)方法。8.3.1龙格—库塔方法的基本思想

我们首先从欧拉公式及改进的欧拉公式着手进行分析。

欧拉公式可改写为

用它计算yn+1需要计算一次f(x，y)的值。若设yn=y(xn)，则yn+1的表达式与y(xn+1)在xn处的泰勒展开式的前两项完全相同，即局部截断误差为O(h2)。改进的欧拉公式又可改写为

用它计算yn+1需要计算两次f(x，y)的值。若设yn=y(xn)，则yn+1的表达式与y(xn+1)在xn处的泰勒展开式的前三项完全相同，即局部截断误差为O(h3)。这两组公式在形式上有一个共同点：都是用f(x，y)在某些点上值的线性组合得出y(xn+1)的近似值yn+1，而且增加计算f(x，y)的次数，可提高截断误差的阶。

因此可考虑用函数f(x，y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xn，yn)处的泰勒展开式与解y(x)在xn处的泰勒展开式的前面几项重合，从而获得达到一定精度的数值计算公式。这就是龙格—库塔的基本思想。8.3.2龙格—库塔方法的推导

按照上述思想，龙格—库塔方法的一般形式设定为

(8.7)

其中，ωi，αi，βij为待定常数。从上式可以看到用它计算yn+1需要计算r次f(x，y)的值，因此(8.7)式称为r阶R-K方法。

下面我们来了解几种常用的R-K方法。

1.二阶龙格—库塔方法

当r=2时，龙格—库塔格式为

适当选取参数ω1，ω2，α2，β21的值，使得在yn=y(xn)的假设下，局部截断误差为y(xn+1)－yn+1=O(h3)。为此把K1，K2代入第一式yn+1=yn+ω1K1+ω2K2中，得然后在(xn，yn)处作泰勒展开，可得

(8.8)

这里记fn=f(xn,yn)。而微分方程y′=f(x,y)的精确解y=y(x)在点xn处的泰勒展开式为

(8.9)把(8.8)式与(8.9)式进行比较，为使y(xn+1)-yn+1=O(h3),令h,h2项系数相等，则有下列方程组:

此方程有无穷多个解，从而有无穷多个二阶龙格—库塔方法。按此方程组解出的每一组解，对应的二阶龙格—库塔方法都是二阶的。常用的二阶龙格—库塔方法有

(1)当ω1=ω2=

,α2=β21=1时，有

这正好是改进的欧拉公式。

(2)当ω1=0,ω2=1,α2=β21=时，有

这种方法称为中间点法。

2.三阶龙格—库塔方法

此时r=3，一般格式为把K1，K2，K3代入yn+1的表达式中，再在(xn，yn)处作泰勒展开，然后与y(xn+1)在xn点处的泰勒展开式作比较，并且要使y(xn+1)－yn+1=O(h4)。可得下列方程组：

此方程有无穷多个解，因此有无穷多个三阶龙格—库塔方法。常用的三阶龙格—库塔方法有

此方法称为三阶库塔方法。

3.四阶龙格—库塔方法

在实际应用中，最常用的龙格—库塔方法是四阶龙格—库塔方法。类似于二、三阶龙格—库塔方法的推导，可得一个含有13个未知量、11个方程的方程组。由于推导复杂，这里从略。此处只介绍最常用的一种四阶龙格—库塔方法：

此方法称为经典四阶龙格—库塔方法。

例8.4取步长h=0.2，用经典四阶R-K公式求解初值问题:

解分析：按照龙格—库塔法的解题步骤应先算出各Ki(i=1，2，3，4)，再代入

中计算。已知f(x，y)=2xy，x0=0，y0=1，h=0.2，由四阶龙格—库塔公式可得

代入yn+1=yn+

(K1+2K2+2K3+K4)(其中n=0,1,2,…，5)求解。本例方程的解为

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

8.4线性多步法

8.4.1线性多步法的基本思想

单步法在计算yn+1时，只用到前一步的信息yn，而没有用到前面几步计算得到的信息。此时要提高精度，需重新计算多个点处的函数值，如R-K方法、计算量较大。若能充分

利用第n+1步前面的多步信息来计算yn+1，便可获得更高精度的算法，这就是多步法的基本思想。多步法中最常用的是线性多步法。如果记y(xn)的近似值为yn，xn=x0+nh，并记fn=f(xn，yn)，则k步线性多步法的一般形式为

(8.10)

其中，αj，βj(j=0，1，…，k)都是实常数，且αk≠0，|α0|+|β0|≠0。一般可设αk=1，此时(8.10)式可写为

(8.11)如果βk≠0，(8.10)式就是隐式方程；如果βk=0，它就是显式方程。如果k=1，(8.10)式就是一种单步法。比如：

k=1,

α0=－1,β0=1,β1=0时，就是欧拉方法；k=1，

α0=－1，β0=，β1=时，就是梯形方法。

下面介绍线性多步法的阶和局部截断误差的定义。

对于一般的线性多步法((8.11)式)，定义算子R为

(8.12)其中,y(x)为区间［a,b］上任意一个连续可微函数。若y(x)充分可微，将y(x+jh)及y′(x+jh)按泰勒式展开有

把它们代入(8.12)式可得

(8.13)其中,(8.14)给出如下定义：

定义8.3若在(8.13)式中

c0=c1=…=cp=0，cp+1≠0

则称算子R［y(x)；h］对应的线性多步法((8.10)式)是p阶

方法。

定义8.4设y(x)是初值问题的解，(8.10)式是一种线性多步法，则

称为是(8.10)式在xn+k处的局部截断误差。Rn+k按h展开的首项称为主局部截断误差。

若(8.10)式是一种p阶线性多步法，则由(8.13)式可知

Rn+k=cp+1hp+1y(p+1)(xn)+O(hp+2)

因此，主局部截断误差就是cp+1hp+1y(p+1)(xn)，而cp+1称为主局部截断误差系数。8.4.2线性多步法的构造

构造线性多步法有多种方法，常用的是数值积分法和泰勒展开法(待定系数法)。

1.数值积分法

初值问题可以写成与之等价的积分方程形式设{xi}为等距节点，xi+1=xi+h。拉格朗日插值基函数li(t)为

(i=0，1，2，…，p)

用拉格朗日插值多项式来近似代替(8.14)式中的被积函数，得到近似计算公式：

若用yn-i代替y(xn－i)，用fn－i表示f(xn－i,yn－i)，则可得线性多步法显式公式:其中,

(i=0,1,2,…，p)

若令t=xn+sh，则

(i=0,1,2,…，p)

当k，j，p取不同值时，得到不同类型的具体公式。如当k=1，j=0，p=0，1，2，…时，可得到阿达姆斯显式

方法：

(i=0，1，2，…，p)

表8.4给出了显式阿达姆斯公式的系数。在阿达姆斯显式方法中，最常用的是p=3的情形：

又如当k=0，j=1，p=0，1，2，…时，得到阿达姆斯隐式公式(用n+1代替n)：

表8.5给出了隐式阿达姆斯公式的系数。在阿达姆斯隐式方法中，最常用的是p=3时的情形：

p=3时，阿达姆斯显、隐式方法都是四阶方法。

2.泰勒展开法(待定系数法)

用泰勒展开法构造线性多步法就是在线性多步法的一般公式

中，假设yn+j=y(xn+j)(j=0，1，2，…，k)，将y(xn+j)和y′(xn+j)在xn处用泰勒公式展开，与推导公式(8.13)式的方法完全类似，导出局部截断误差按h的升幂排列表达式，然后

再确定相应的待定系数αj，βj。

我们举几种特殊情形给予说明。考虑线性两步法，此时k=2，α2=1，则有

yn+2+α1yn+1+α0yn=h(β2fn+2+β1fn+1+β0fn)

由(8.14)式可得系数α0,α1,β0,β1,β2满足:解此方程可得

从而一般线性二步方法可写为

(8.15)对任意α0，有当α0≠－1时，c4≠0，(8.15)式是三阶的；当α0=－1时，c4=0,c5≠0，(8.15)式就是二步四阶公式(辛普森公式)：

其截断误差为又如当α0=－5时，方法是显式的，否则就称为隐式的，此时(8.15)式称为二步四阶米尔恩方法。

考虑线性三步法，此时k=3，α3=1，则有由(8.14)式可得系数α0,α1,α2,β0,β1，β2，β3满足的方程为

此方程组中含有六个未知量、四个方程，因此有无穷多组解，也即有无穷多个三步方法。如令α0=α1=0，可得到

此时可得到阿达姆斯三阶显式方法:

其局部截断误差为又如,令α0=α1=0,α2=－1，且考虑隐式方法，可得到阿达姆斯四阶隐式方法:

其截断误差为类似地，还可得到其他线性多步方法，如米尔恩方法:

哈明(Hamming)方法:

在应用线性多步法求解初值问题时，开始几点处的函数值往往用别的方法计算。一般地选用与多步法同阶的单步法，如龙格—库塔方法、泰勒方法等。对线性隐式多步法，除开始几点的函数值需单独计算外，还需迭代求解或采用预估—校正方法求解。

例8.5应用四步四阶阿达姆斯显式方法求解初值问题：

取步长h=0.1。

解分析：应用四步显式法必须有四个起步值，y0已知，而y1,y2，y3可用精度相同的四阶龙格—库塔方法求出。

步长h=0.1，节点xn=nh=0.1n(n=0,1,…，6)，由四阶龙格—库塔法算出：

y1=1.0048375,

y2=1.0187309,

y3=1.0408184

四步四阶阿达姆斯显式为已知fn=f(xn,yn)=xn－yn+1，h=0.1,xn=0.1n，则

由此算出

y4=1.0703231,

y5=1.1065356,

y6=1.1488186

例8.6以四步四阶阿达姆斯显式与隐式作为预估—校正公式解初值问题：

取步长h=0.1，并将结果与精确解比较。

解分析：由于预估公式与校正公式都是四阶的，因此起步值也应按四阶公式求出。

已知y0=1，按四阶龙格—库塔公式法算出：

y1=1.095446，y2=1.183217,

y3=1.264912

四阶阿达姆斯显式作为预估式(取步长h=0.1)，即四阶阿达姆斯隐式作为校正式，即

该问题的精确解为。

此题的数值计算结果和精确结果及误差列于表8.6。

8.5算法的稳定性及收敛性

8.5.1算法的稳定性

稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法求解过程中存在着各种计算误差，这些计算误差，如舍入误差等引起的扰动在传播

过程中，可能会大量积累，对计算结果的准确性将产生影响，这就涉及到算法稳定性问题。本节中仅介绍绝对稳定性的概念。

定义8.5当在某节点xn上的yn值有大小为δ的扰动时，如果在其后的各节点xj(j>n)上的值yj产生的偏差都不大于δ，则称这种方法是绝对稳定的。

稳定性不仅与算法有关，而且与方程中函数f(x,y)也有关，讨论起来比较复杂。为简单起见，通常只针对模型方程

y′=λy

(λ<0)

(8.16)

来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程如此简单的问题都不稳定的话，也就很难用于其他方程的求解。下面以欧拉方法为例讨论绝对稳定性。

先考察向前欧拉方法的稳定性。模型方程

y′=λy

(λ<0)

的欧拉公式为

yn+1=yn+λhyn=(1+λh)yn

设yn有误差δn，则实际参与运算的，由此引起yn+1的误差为δn+1，实际得到

即有

yn+1+δn+1=(1+λh)(yn+δn)从而有

δn+1=(1+λh)δn

要使|δn+1|<|δn|，必须有|1+λh|<1。因此向前欧拉方法的绝对稳定域为|1+λh|<1。在复平面上，|1+λh|<1是以1为半径，以－1为圆心的圆内部，所以向前欧拉方法的绝对稳

定域是圆域，如图8-2所示。图8-2向前欧拉方法的绝对稳定域用向后欧拉方法解模型方程(8.16)式的计算公式为

yn+1=yn+λhyn+1

解出，则对于向后欧拉公式法，δn应满足。由于λ<0，则有，故恒有|δn+1|<|δn|，因此向后欧拉方法是绝对稳定的，它的绝对收敛域为|1－λh|>1。在复平面上，它是以1为中心，以1为半径的圆外部，如图8-3所示。图8-3由图8-2和图8-3可知，向后欧拉方法的绝对稳定域比向前欧拉方法的绝对稳定域大得多，但可以证明这两种方法的收敛阶数是相同的，只是向前欧拉方法是显式方法，向后欧拉方法是隐式方法。这也说明，隐式方法的稳定性一般比同阶的显式方法的稳定性要好得多。

用二级二阶龙格—库塔方法求解试验方程(8.16)式的计算公式为

yn+1=yn+h［ω1λyn+ω2λ(yn+β21hyn)］

=［1+(ω1+ω2)λh+ω2β21(λh)2］yn利用二级二阶龙格—库塔方法的参数可知：

类似于向前欧拉方法的分析，可得二级二阶龙格—库塔方法的绝对稳定域为

也就是|1+λh|<1，它与向前欧拉方法的绝对稳定域相同。同理可得三级三阶龙格—库塔方法的绝对稳定域是

绝对稳定区间是(－2.51,0)。四级四阶龙格—库塔方法的绝对稳定域是

绝对稳定区间是(－2.78,0)。

二级二阶、三级三阶及四级四阶龙格—库塔方法的绝对稳定域如图8-4所示。图8-48.5.2算法的收敛性

微分方程初值问题数值解法的基本思想是将微分方程转化为差分方程来求解，并用计算值yn近似代替y(xn)。这种近似代替是否合理，还须看分割区间［xi－1，xi］的长度h越来越小时，即h=xi－xi－1→0时，yn→y(xn)是否成立。若成立，则称该方法是收敛的；否则称为不收敛。这里仍以欧拉方法为例，来分析收敛性。

欧拉格式如下：

yn+1=yn+hf(xn，yn)

设表示取yn=y(xn)时，按欧拉公式的计算结果，即

欧拉方法局部截断误差为设有常数，则

(8.17)

总体截断误差为

又

(8.18)由于f(x,y)关于y满足李普希兹条件，即有

|f(xn,y(xn))－f(xn,yn)|≤L|y(xn)－yn|

(8.19)

将(8.19)式代入(8.18)式，有

再利用(8.17)式、(8.18)式可得到

(8.20)即

|εn|≤(1+hL)|εn－1|+ch2

上式反复递推后，可得设xn－x0=nh≤T(T为常数)，因为1+hL≤ehL，故(1+hL)n

≤enhL≤eTL把上式代入(8.20)式，得

若不计初值误差，即ε0=0，则有

(8.21)

(8.21)式说明，当h→0时，εn→0，从而yn→y(xn)，所以欧拉方法是收敛的，且其收敛速度为O(h)，即具有一阶收敛速度。

8.6一阶常微分方程组与高阶方程

前面已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。8.6.1一阶常微分方程组

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

可以把单个方程y′=f(x,y)中的f,y看做向量来处理，这样就可把前面介绍的各种算法推广到求一阶方程组初值问题中来。设xn=x0+nh(n=1,2,3,…)，yn，zn为节点xn上的近似解，则有下列求解格式：

(1)向前欧拉格式：

这是一步一阶方法，为显式。

(2)改进的欧拉格式：

这是一步二阶方法。

(3)四阶标准龙格—库塔格式：

(8.22)其中，

(8.23)把节点xn上的yn和zn值代入(8.23)式，依次求出K1,L1,K2,L2,K3,L3,K4,L4，然后把它们代入(8.22)式，算出节点xn+1上的yn+1和zn+1值。

对于具有三个或三个以上方程的方程组的初值问题，也可用类似方法处理。此外，多步法也同样可应用于求解方程组初值问题。

例8.7用改进的欧拉方法求解初值问题：

取步长h=0.1，保留6位小数。

解改进的欧拉公式为将f(xn,yn,zn)=xy－z,g(xn,yn,zn)=及h=0.1代入上式，得由初值y0=y(0)=1,z0=z(0)=2计算得8.6.2高阶微分方程

高阶微分方程(或方程组)的初值问题可以通过变量代换化为一阶方程组来求解。例如，有二阶微分方程的初值问题:

(8.24)作变换z=y′，则(8.24)式可化为一阶方程组的初值问题:

(8.25)

对此问题就可用上小节的方法求解。此方法同样可以用来处理三阶或更高阶的微分方程(或方程组)的初值问题。

例8.8采用改进的欧拉方法求解二阶常微分方程初值问题:

(0≤x≤1)

取步长h=0.1，计算y(0.1)的近似值(最后结果保留小数点后

5位)。

解令z=y′，则所给微分方程初值问题等价于解该一阶微分方程组，改进的欧拉公式为代入h=0.1,y0=－0.4,z0=－0.6,x0=0.0，计算得

由此可得

y(0.1)≈y1=－0.46200

例8.9应用四阶龙格—库塔方法求解二阶常微分方程初值问题:

取步长h=0.1，计算y(0.1),z(0.1)的近似值(最后结果保留小数点后5位)。

解令z=y′，则所给微分方程初值问题等价于解该一阶微分方程组，四阶龙格—库塔公式为其中,取步长h=0.1,x0=0,y0=0,z0=1,代入上式得则

由此可得

y(0.1)≈y1=0.10369,

z(0.1)≈z1=1.11043

8.7解微分方程的波形松弛方法

本节介绍波形松弛方法(也叫动力迭代方法)在求解微分方程的初值问题和边值问题方面的应用。

通过积分方法可以求出一些特殊类型微分方程的解析解，但对维数较高而且形式复杂的微分方程，或者无法求得其解析解，或者所得到的解析解形式复杂难以得到所求解的性质的问题，波形松弛方法是求解它们的一个有效工具。求解微分方程初值问题及边值问题的迭代方法较多，本节仅介绍在毕卡(Picard)逐次逼近的基础上发展起来的波形松弛方法在求解微分方程初值问题以及边值问题方面的简单应用。8.7.1微分方程初值问题的波形松弛方法

考虑微分方程的初值问题：

(8.26)

其中，f：R×Rn→Rn，y∈Rn,f满足解的存在唯一性条件。为了叙述清楚，我们先以二阶自治系统情形为例。

(8.27)

类似于线性代数方程组迭代解法中的迭代格式，有如下动力迭代格式。

1)雅可比动力迭代

雅可比动力迭代的基本迭代格式为

(k=0，1，2，…；x∈［x0,T］)

(8.28)通过逐次迭代，产生一系列函数序列{(y1(k)(x)，y2(k)(x))T},如果

则(y1(x)，y2(x))T就是(8.27)式的解。

雅可比动力迭代的优点在于：对每一个k值，如果第k步已求出y1(k)(x)和y2(k)(x)，那么在进行第k+1步的求解过程中，(8.28)式中每个方程可以分别独立地求解。这种分解格式特别对高维微分方程具有很大的优势，因为它可以进行并行计算，只在求出结果时，互相传递信息。

2)高斯—赛德尔动力迭代

高斯—赛德尔动力迭代的格式为

(k=0，1，2，…；x∈［x0,T］)

(8.29)同样，经过逐步迭代可得解序列{(y1(k)(x)，y2(k)(x))T}，当

时，就可得到(y1(x)，y2(x))T，它是(8.28)式的解。高斯—赛德尔动力迭代的优点在于由(8.29)式中第一个方程求出解y1(k+1)(x)时，求解第二个方程时就用y1(k+1)(x)，而不用y1(k)(x)。这是因为在{(y1(k)(x),y2(k)(x))T}收敛时，y1(k+1)(x)比y1(k)(x)更接近于y1(x)，因此，在收敛情况下，我们把第一个方程已经算出的y1(k+1)(x)代入第二个方程,同时也能使第一个方程的解序列收敛更快。高斯—赛德尔动力迭代的缺点在于只有当第一个方程求出y1(k+1)(x)时，才能计算第二个方程，因此不能在并行机上进行计算。例如,应用高斯—赛德尔动力迭代方法来计算下面两个变量的微分方程初值问题的解：

初始函数选为，迭代一次有迭代两次有

迭代三次有迭代n次有显然，迭代序列{(y1(k)(t)，y2(k)(t))T}收敛到{(sint，cost)T},

即

y1(∞)(t)=sin(t)，y2(∞)(t)=cos(t)

事实上，我们容易验证该二阶微分方程的解为

y1(t)=sin(t)，y2(t)=cos(t)

迭代序列收敛到精确解的过程如图8-5和图8-6所示，图中给出了n=0，1，2，3，4，5时y1(n)(t)，y2(n)(t)的曲线图。图8-5图8-6

3)SOR动力迭代

SOR动力迭代的基本格式为

(8.30)如果迭代函数序列收敛，即

则称(y1(x)，y2(x))T为(8.27)式的解。

SOR迭代过程是把前一次得到的结果记为(z1(k)(x),z2(k)(x))T,然后把两次所得结果进行加权平均，权系数为ω，所得结果作为下一次迭代的初值。SOR动力迭代的优点是可以使迭代序列收敛速度更快，缺点是计算量大，同时权系数ω(也称为松弛因子)的选取比较困难。对一般形式的n维非自治微分动力系统((8.26)式)，其动力迭代格式为

(8.31)此处，G和g是一般的分裂函数，且满足：

G：［x0，T］×Rn×Rn×Rn→Rn，

G(x，y，y，y)=f(x，y)

g：［x0，T］×Rn×Rn→Rn，

g(x，y，y)=f(x，y)

事实上，(8.31)式的迭代格式包括了微分方程中毕卡逐次逼近的思想，它的一个更简单的形式为

(8.32)此处，函数F满足：

F：［x0，T］×Rn×Rn→Rn，

F(x，y，y)=f(x，y)

我们把(8.26)式中函数f的n个分量记为f1，f2，…，fn，且n维向量y(k)的n个分量记为y1(k)，y2(k)，…，yn(k)，则(8.26)式的雅可比、高斯—赛德尔、SOR动力迭代格式按分量写出来分别为

(8.33)

(8.34)

(8.35)应当注意，在给出上述迭代格式以及迭代函数序列的收敛时没有考虑收敛对自变量取值范围的一致性问题。在应用动力迭代方法求解微分系统问题时一般要注意以下三个问题：

(1)积分区间上的误差不是一致的。初始函数在起始阶段能有比较好的逼近效果，但是随着自变量x的增加，情况就会不断恶化。

(2)每步迭代延长波形逼近真实解所需要的时间。

(3)计算波形和真实解之间的误差未必在每个时间节点减少。从上一步迭代到下一步迭代，计算误差反倒可能增加，尤其对于时间值x较大的范围。8.7.2微分方程初值问题波形松弛方法的收敛问题

下面讨论线性微分方程初值问题动力迭代的收敛性。收敛性结论主要基于压缩映像原理。

线性微分方程初值问题为

(8.36)其中，A=(aij)n×n为已知矩阵；g(x)为n维向量函数，

g：R→Rn，y∈Rn,x∈R。对初值问题(8.36)式，其解析解为

而(8.36)式的动力迭代格式为

(8.37)其中，A=M－N。此处，M为对角矩阵或块对角矩阵；N为非对角耦合矩阵。求解(8.37)式可得：

(8.38)

令

定义算子K为则(8.38)式可写成不动点迭代形式：

y(k+1)(x)=Ky(k)(x)+Φ(x)(8.39)

容易证明(8.39)式的不动点是(8.36)式的解，且(8.36)式的解是(8.39)式的不动点。以y(x)表示(8.36)式的解，令ε(k)(x)=y(k)(x)－y(x)，则由(8.39)式可得

ε(k+1)(x)=Kε(k)(x)

从而有

ε(k)(x)=Kkε(0)(x)

(8.40)

此处Kk表示算子K的k重卷积。显然,(8.40)式收敛的充分必要条件为ρ(K)<1。为了获得第k步误差ε(k)(x)对初始误差ε(0)(x)的真实估计，首先给出如下定义。

定义8.6令‖·‖表示Cn中的任一种确定范数，则在［0,T］上，连续函数的最大范数定义为由定义8.6可以证明：

(8.41)

此处,表示算子γ的k重卷积γ*γ*…*γ，且γ*u(x)=

。由于已给T<+∞，因此有

‖γ‖T≤c

(8.42)且c=eT‖M‖‖N‖，则

由归纳可得

结合(8.39)式就有定理8.3。

定理8.3由动力迭代算法((8.37)式)以及K和γ(x)的定义，并由(8.42)式可得算子K的一个界为

(8.43)

注：由于，因此定理8.3表明对所有有限T值，由(8.37)式或(8.39)式产生的函数序列{y(k)(x)}收敛，收敛的速度可由K的分解表达式进行分析。由于

因此对任何固定T<+∞和非零λ，Rλ(K)是有界的。8.7.3微分方程边值问题的波形松弛方法

微分方程边值问题的基本形式为

(8.44)

对边值问题的动力迭代，其基本格式完全类似于对初值问题的迭代。只要把迭代格式(8.33)式、(8.34)式和(8.35)式中相应的初值条件y(k+1)(x0)=y0改为y(k+1)(0)=y(k+1)(T)即可得到微分方程边值问题的雅可比、高斯—赛德尔、SOR动力迭代格式。我们在此仅介绍微分方程边值问题的动力迭代方法中的一个收敛性条件。

考虑微分方程边值问题(8.44)式的动力迭代格式：

(x∈［0，T］；k=0，1，2，…)

(8.45)

其中，F：［0，T］×Rn×Rn→Rn，满足F(x,y,y)=f(x,y)，

且x(0)(t)是一个初始猜测函数，满足x(0)(0)=x(0)(T)。令‖·‖表示Rn中的2—范数，〈·，·〉表示标准内积，即〈w，w〉=‖w‖2。假设F满足下列条件：

(1)存在非负函数L1(x)，使得

‖F(x,u1,v)－F(x,u2,v)‖≤L1(x)‖u1－u2‖

(8.46)

其中，x∈［0,T］,v∈Rn,u1,u2∈Rn

。

(2)存在函数a(x)，使得

〈F(x,u,v1),F(x,u,v2)〉≤a(x)‖v1－v2‖

(8.47)

其中,x∈［0,T］,u∈Rn,v1,v2∈Rn,且。

先给出引理8.1。

引理8.1设a(x)∈R,Q(t)∈R,x∈［0,T］,如果φ(x)

满足:

则有

证明首先证明当Q(x)≡0,x∈［0,T］时，φ0=0结论成立，因为

所以有即φ(x)≤0，x∈［0，T］。其次我们来证明其一般情形。

令

(x∈［0，T］)

则φ(x)满足下列初值问题:我们定义

则有即

由前述情形可知，y(x)≤0，x∈［0，T］，即φ(x)≤ψ(x)，x∈［0，T］，从而完成了证明。

定理8.4对(8.44)式和(8.45)式，有

(8.48)

证明首先，令ε(k+1)(x)=y(k+1)(x)－y(x)，由(8.46)式和(8.47)式，有又由于

则对‖ε(k+1)(x)‖≠0,我们有

(8.49)

由引理8.1可知：

(8.50)此不等式对‖ε(k+1)(x)‖=0也成立。由于‖ε(k+1)‖(0)=

‖ε(k+1)‖(T)，因此由(8.50)式有

(8.51)

把(8.51)式代入(8.50)式，则有(8.48)式成立。

由引理8.1和定理8.4，可以得到定理8.5。

定理8.5对微分方程边值问题(8.44)式和(8.45)式，如果函数f满足条件(8.46)式和(8.47)式，且

则动力迭代序列{y(k)(x)}收敛到(8.44)式的解。

证明由于

定义范数‖·‖e如下：此处，‖·‖是Rn

中的范数，且λ是正常数，而

是有界的，因为此表达式内的函数在［0，T］上是连续的，且

(8.52)

此处，因此我们有

(8.53)

此处

根据条件(8.52)式及定理假设条件，可以选取充分大λ，使L<1，即就是迭代格式(8.53)式收敛，从而{ε(k+1)(x)}收敛到零，证毕。下面我们就L1(x)=L1为常数，且a(x)=L2<0，且为常数给出一个数值计算的例子。

例8.10考虑如下二维系统：

此处，,我们取如下动力迭代格式:

(x∈［0,2π］；k=0,1,…)

表8.7给出了迭代次数与误差的数值。8.8微分方程边值问题的数值方法

由于高阶微分方程经过变量代换可以化为一阶微分方程组，因此微分方程组边值问题的一般形式为

其中，F：R×Rn→Rn，y∈Rn。假设F(x，y)满足解的存在唯一性条件。在实际工程技术中经常遇到二阶微分方程：

y″=f(x,y,y′)

(x∈［a,b］)

(8.54)

为了确定唯一解，需要附加两个定解条件。当定解条件取为解在区间［a,b］两端点的状态时，相应的定解条件称为两点边值问题。

边值条件有以下三类提法。

第一类边界条件：

y(a)=α,

y(b)=β

当α=0或β=0时，称为齐次边界条件；否则称为非齐次边界条件。第二类边界条件：

y′(a)=α,

y′(b)=β

当α=0或β=0时，称为齐次边界条件；否则称为非齐次边界条件。

第三类边界条件：

y(a)－α0y′(a)=α1,

y(b)－β0y′(b)=β1

其中，α0≥0，β0≥0，α0+β0>0。当α1=0或β1=0时，称为齐次边界条件；否则称为非齐次边界条件。微分方程(8.54)式附加第一、二、三类边界条件后，分别称为第一、二、三类边值问题。

下面仅介绍二阶微分方程两点边值问题的打靶法和有限差分方法。8.8.1打靶方法

以二阶微分方程第一类边值问题为例来讨论打靶法。它的基本原理是将边值问题转化为相应的初值问题：

(8.55)

来求解。令y1=y′，可以把(8.55)式化为

(8.56)至此，问题转化为求适当的z，使初值问题(8.56)式的解y(x,z)在x=b处的值满足右端边界条件

y(b,z)=β

这样,初值问题(8.56)式的解y(b,z)就是相应边值问题的解。而求初值问题(8.56)式可用前面介绍的数值方法来求解，比如牛顿方法或其他迭代方法。先用线性插值方法在z=z0,

即y′(a)=z0

时求解初值问题(8.55)式，得解y(b,z0)=β0。若|β－β0|≤ε(ε为允许误差)，则y(xj，z0)(j=0，1，2，…，m)是初值问题(8.56)式的数值解，也就是相应边值问题的解。若|β－β0|>ε，则调整初始条件为y′(a)=z1，重新解初值问题(8.55)式，得解y(b,z1)=β1。若|β－β1|≤ε(ε为允许误差)，则y(xj,z1)(j=0,1,2,…，m)为所求，否则再修正z1为z2