第七章 常微分方程模型的数值解法

第七章 常微分方程数值解法简介微分方程在科学和工程技术中有很广泛的应用。

许多实际问题的数学模型都可以用微分方程来描述，归结为常微分方程的定解问题；很多偏微分方程问题，也可以化为常微分方程问题来近似求解，但是求出所需的解绝非易事。

实际上，除了极特殊情形外，人们不可能求出微分方程的解析解，只能用各种近似方法得到满足一定精度的近似解。

在常微分方程中已经熟悉了级数解法和Picard 逐步逼近法，这些方法可以给出解的近似表达式，称为近似解析方法。

另一类方法只给出解在一些离散点上的值，称为数值方法。

后一类方法应用范围更广，特别适合用计算机计算，本章主要介绍常用的常微分方程数值解法。

7.1实际问题的微分方程模型函数是事物的内部联系在数量方面的反映,如何寻找变量之间的函数关系,在实际应用中具有重要意义。

在许多实际问题中，往往不能直接找出变量之间的函数关系,但是有时却容易找出变量的改变量之间的关系，从而建立描述问题的微分方程模型。

例7.1.1 将初始温度00150u C =的一碗汤放置于环境温度a u 保持为024C 的桌上，10分钟后测得汤的温度为0100C 。

如果汤的温度低于055C 才可以喝，试问再过20分钟后这碗汤能喝了吗？解：为了解决这一问题，需要了解有关热力学的一些基本规律。

热量总是从温度高的物体向温度低的物体传导的；在一定的温度范围内，一个物体的温度变化速度与这个物体的温度和其所在的介质温度的差值成正比。

设物体在t 时刻的温度为()u u t =，从t t t →+∆温度从()()u t u t t →+∆，注意到热量总是从温度高的物体向温度低的物体传导，因而0a u u >，所以温度差a u u -恒正，又因物体将随时间而逐渐冷却；则温度的改变量为:()()(())a u u t t u t k u t t u t∆=+∆-=-+∆-∆两边除以t ∆，并令0t ∆→得温度变化速度为：()a du k u u dt=--这里0k >是比例常数。

从而得出描述物体冷却过程的微分方程模型为：0()(0)a duk u u dt u u ⎧=--⎪⎨⎪=⎩（7.1.1）容易求出这个一阶微分方程初值问题的解为：0()kta a u u u u e-=+- （7.1.2）根据所给问题的条件知当10t =时，01100u u C ==，得到：1010()ka a u u u u e-=+-将00150u C =，024a u C =带入得到：0111ln ln 1.660.0511010a au u k u u -=≈≈-从而得这碗汤的温度随时间变化的函数关系为： 0.0510.0510()()24126tta a u u t u u u ee--==+-=+ （7.1.3）于是，将30t =带入计算得，再过20分钟后汤的温度就是0251u C ≈。

说明再过20分钟后这碗汤能喝了。

不过并不是所有的微分方程模型都可求出解析解，从而得出实际问题的解释。

例如看似简单的微分方程22dy x ydx=+，自德国数学家Wilhelm von Leibniz 提出100多年后才被法国数学家Joseph Liouville 证明它没有解析解，只能借助于数值的方法求数值解。

例7.1.2 某地区发现一种有免疫性的传染病，为了控制疫情扩散对该地区人群进行隔离处理。

为了分析受感染人数的变化规律，需要建立描述传染病传播过程的数学模型。

解：设该地区的总人数为常数N ，任意t 时刻病人、健康人和病人治愈后移出感染系统的移出者的比例分别为(),(),()i t s t r t ，病人的日接触率λ，日治愈率μ，则容易得出从t t t →+∆时刻，病人和健康人的改变量为：[()()]()()()[()()]()()N i t t i t N s t i t t N i t tN s t t s t N s t i t t λμλ+∆-=∆-∆⎧⎨+∆-=-∆⎩每个方程两边除以t ∆，并令0t ∆→，化简后得：00(0),(0)d i s i i d td s s i d t i i s s λμλ⎧=-⎪⎪⎪=-⎨⎪==⎪⎪⎩（7.1.4） 其中， ()()()1s t i t r t ++=，对任意的t 。

（7.1.4）就是描述病人和健康人的比例()i t 和()s t 随时间变化的微分方程模型，这是一个微分方程组的初值问题。

但是这一初值问题的解析解是无法求出的，因此不能直接利用()i t 和()s t 的解析式来分析和解决问题。

在数学建模课程中学到的大量数学模型都是微分方程形式给出的，各类微分方程本身和它们的解所具有的特性已在常微分方程及数学物理方程中得以解释。

虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程，在实际应用中人们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值。

这样一组近似解称为微分方程在该范围内的数值解，寻找微分方程数值解的过程称为微分方程的数值解法。

7.2简单的数值方法与基本概念7.2.1常微分方程初值问题设(,)f x y 在区域:,G a x b y ≤≤<∞上连续，求()y y x =满足0(,),()dyf x y a x b dx y a y⎧=<≤⎪⎨⎪=⎩（7.2.1） 其中0y 是已知常数,这就是一阶常微分方程的初值问题。

为使问题（7.2.1）的解存在、唯一且连续依赖初值0y ，即初值问题（7.2.1）适定，还必须对右端项(,)f x y 加适当限制，通常要求(,)f x y 关于y 满足是已知函数,且满足Lipschitz 条件：存在常数L ，使1212(,)(,)f x y f x y L y y -≤- （7.2.2）对所有[,]x a b ∈和12,(,)y y ∈-∞+∞成立。

本章总假定f 满足此条件。

7.2.2 Euler 法及改进的Euler 法 7.2.2.1 Euler 方法的导出与几何意义最简单的数值解法是Euler 法。

将区间[0,]T 作N 等分，小区间的长度/h T N =称为步长，点列i x a ih =+，(0,1,...,)i N =称为节点，0x a =。

由已知初值00()y x y =，可算出()y x 在0x x =的导数'00000()(,())(,)y x f x y x f x y ==。

下面用三种方法导出Euler 法。

本章用()i y x 表示函数()y x 在i x 点的精确值、i y 表示()i y x 的近似值。

（1）幂级数展开法利用Taylor 展式23''''''100000000()()()()()()26(,)hhy x y x h y x hy x y x y y hf x y R ξ=+=+++=++ （7.2.3）其中01(,)x x ξ∈，并略去二阶小量0R ，得1000(,)y y hf x y =+，1y 就是1()y x 的近似值。

利用1y 又可算出2y ，如此下去可算出()y x 在所有节点上的值，一般递推公式为：1(,)n n n n y y hf x y +=+，0,1,...,1n N =- （7.2.4）这就是Euler 法。

Euler 法有明显的几何意义。

实际上，初值问题（7.2.1）的解是x ，y 平面上过点00(,)x y 的一条积分曲线。

按Euler 法，过初始点00(,)x y 作经过此点的积分曲线的切线（斜率为00(,)f x y ），沿切线取点11(,)x y （1y 按（7.2.4）计算）作为11(,())x y x 的近似。

然后过11(,)x y 作经过此点的积分曲线的切线（斜率为11(,)f x y ），沿切线取点22(,)x y （2y 按（7.2.4）计算）作为22(,())x y x 的近似。

如此下去，即得一以(,)n n x y 为顶点的折线，这就是用Euler 法得到的近似积分曲线（如图7.1.1所示）。

从几何上看，h 越小，此折线逼近积分曲线越好，因此也称Euler 法为Euler 折线法。

图7.1.1 Euler 近似积分曲线（2）数值微分法 利用向前差商近似导数1()()()n n n y x y x y x h+-'≈1()()()(,)n n n n n n y x y x hy x y h f x y +'≈+=+从而得出Euler 法的一般递推公式：1(,)n n n n y y hf x y +=+，0,1,...,1n N =-（3）数值积分法将初值问题（7.2.1）写成等价的积分形式：0()()(,())x x y x y x f t y t dt =+⎰取1x x =时得1010()()(,())x x y x y x f t y t dt =+⎰（7.2.5）用左矩形公式近似右端积分，并用1y 替代1()y x ，即得1000(,)y y hf x y =+，从而也可得出Euler 法的一般递推公式为：1(,)n n n n y y hf x y +=+，0,1,...,1n N =-。

7.2.2.2 改进的Euler 方法由Euler 方法的数值积分导出法可知只要给出（7.2.5）式右端定积分的一种近似计算方法，就可得出初值问题（7.2.1）的一种数值求解方法。

如果用右矩形公式近似（7.2.5）式右端积分，则可得1011(,)y y hf x y =+，从而也可得出一般递推公式为：111(,)n n n n y y hf x y +++=+，0,1,...,1n N =- （7.2.6）称（7.2.6）为后退Euler 法。

如果用梯形公式近似（7.2.5）式右端定积分，则可得100011[(,)(,)]2h y y f x y f x y =++从而得出一般递推公式为：111[(,)(,)]2n n n n n n h y y f x y f x y +++=++，0,1,...,1n N =- （7.2.7）称（7.2.7）为改进的Euler 法，显然改进的Euler 法比Euler 法精度更高。

后退Euler 法和改进的Euler 法，由于未知数1n y +同时出现在等式的两边，故称为隐式；如果未知数1n y +由已知量直接计算即不只出现在等式右端，则称为显式。

对于隐式算法每步计算需要解关于1n y +的方程，而这样的方程往往是非线性的，通常将初值取为[0]1n n y y +=用迭代法求解，一般只需迭代几步即可收敛。

一般先用显式公式计算一个初值，再用隐式公式迭代求解。

如果先用显式Euler 公式作预测，算出1(,)n n n n y y hf x y +=+，再将1n y +代入隐式梯形公式的右边作校正，得到111[(,)(,)]2n n n n n n h y y f x y f x y +++=++。