第 7 章：常微分方程的数值解法（一）

主要内容

• 7.1 引言
• 7.2 数值方法的基本思想
• 7.3 欧拉法
• 7.4 龙格-库塔法
• 7.5 亚当姆斯法
• 7.6 算法的稳定性及收敛性
• 7.7 一阶常微分方程与高阶方程

7.1 引言

包含自变量、未知函数及未知数的导数或微分的方程称为微分方程。在微分方程中，自变量的个数只有一个的，称为常微分方程；自变量的个数有两个或以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数\(y\)及其各阶导数

\[y',y'',...,y^{(n)} \]

都是一次的，则称它是线性的，否则称为非线性的

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

从实际问题中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程的初值问题

\[\left\{ \begin{array}{l} y' = f(x, y) \\ y(x_0) = y_0 \end{array} \right.\]

7.2 数值方法的基本思想

\[\left\{ \begin{array}{l} y' = f(x, y) \\ y(x_0) = y_0 \end{array} \right.\]

• 对如上常微分方程初值问题的数值解法，就是要算出精确解\(y(x)\)在区间\([a,b]\)上的一系列离散节点\(a=x_0<x_1<...<x_{(n-1)}<x_n=b\)处的函数值\(y(x_0),y(x_1),...,y(x_n)\)的近似值\(y_0,y_1,...,y_n\)。相邻的两个节点间距\(h=x{i+1}-x_i\)称为步长，步长可以相等，也可以不等。本章假定步长\(h\)为定值，称为定步长，这时的节点可表示为

\[x_i=x_0+ih,i=1,2,...,n \]

• 数值解法需要把连续性问题加以离散化，从而求出离散节点的数值解
• 对常微分方程数值解法的基本出发点是离散化。其数值解有一个基本特点，它们都采用“先进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息\(y_i,y_{i-1},y_{i-2},...,y_0\)计算\(y_{i+1}\)的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题

\[\left\{ \begin{array}{l} y' = f(x, y) \\ y(x_0) = y_0 \end{array} \right.\]

• 中的导数\(y'\)进行不同的离散化处理
• 单步法：计算\(y_{i+1}\)时只用到\(x_{i+1},x_i,y_i\)，即前一步的值,代表是龙格-库塔法
• 多步法：计算\(y_{i+1}\)用到前几步的值，代表是亚当斯法

7.3 欧拉法

欧拉公式

欧拉方法是解初值问题的最简单的数值方法。初值问题

\[\left\{ \begin{array}{l} y' = f(x, y) \\ y(x_0) = y_0 \end{array} \right.\]

的解\(y=y(x)\)代表通过点\((x_0,y_0)\)的一条称之为微分方程的积分曲线。积分曲线上每一点\((x,y)\)的切线的斜率\(y'(x)\)等于函数\(f(x,y)\)在这点的值

欧拉法的求解过程是：从初始点\(P_0\)（即点\((x_0,y_0)\)）出发，作积分曲线\(y=y(x)\)在\(p_0\)点上的切线\(\overline{P_0P_1}\)（其斜率为\(y'(x_0)=f(x_0,y_0)\)），与直线\(x=x_1\)相交于\(p_1\)点（即点\((x_1,y_1)\)），得到\(y_1\)作为\(y(x_1)\)的近似值，如图所示

• 过点(x_0,y_0)，以\(f(x_0,y_0)\)为斜率的切线方程为

\[y=y_0+f(x_0,y_0)(x-x_0) \]

• 当\(x=x_1\)时，得\(y_1=y_0+f(x_0,y_0)(x_1-x_0)\)
• 这样就得到了\(P1\)点的坐标
• 这样，从\(x_0\)逐个算出\(x_1,x_2,...,x_n\)对应的数值解\(y_1,y_2,...,y_n\)
• 从图形上看，就获得了一条近似于曲线\(y=y(x)\)的折线\(\overline{P_1P_2P_3...P_n}\)

欧拉法的几何意义

通常取\(x_{i+1}-x{i}=h_i=h\)，则欧拉法的计算格式

\[\begin{cases} y_{i+1} = y_i + hf(x_i, y_i) \\ y_0 = y(x_0) \end{cases} \quad i = 0, 1, \ldots, n\]

还可以用数值微分、数值积分和泰勒展开法推导欧拉格式，

以数值积分为例进行推导
将方程\(y'=f(x,y)\)的两端在区间\([x_i,x_{i+1}]\)上积分得

\[\int_{x_{i}}^{x_{i+1}} y' \, dx = \int_{x_{i}}^{x_{i+1}} f(x, y) \, dx \]

\[y(x_{i+1}) = y(x_i) + \int_{x_i}^{x_{i+1}} f(x, y) \, dx = y(x_i) + \int_{x_i}^{x_{i+1}} f[x, y(x)] \, dx \]

选择不同的计算方法计算上式的积分项\(\int_{x_i}^{x_{i+1}} f[x, y(x)] \, dx\)，就会得到不同的积分公式

• 用左矩形方法计算积分项

\[\int_{x_i}^{x_{i+1}} f[x, y(x)] \, dx \approx (x_{i+1} - x_i) f[x_i, y(x_i)] \]

• 带入先前的式子，并用\(y_i\)近似代替式中\(y(x_i)\)即可得到向前欧拉公式

\[y_{i+1} = y_i + h f(x_i, y_i) \]

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

例题1

用欧拉法解初值问题

\[\begin{cases} y' = -y - xy^2 & (0 \leq x \leq 0.6) \\ y(0) = 1 \end{cases}\]

取步长\(h=0.2\)，计算过程保留4位小数

梯形公式

为了提高精度，对方程\(y'=f(x,y)\)的两端在区间\([x_i,x_i+1]\)上进行积分得到

\[y(x_{i+1}) = y(x_i) + \int_{x_i}^{x_{i+1}} f[x, y(x)] \, dx \]

改用梯形方法计算其积分项，即

\[\int_{x_i}^{x_{i+1}} f[x, y(x)] \, dx \approx \frac{x_{i+1} - x_i}{2} \left[ f(x_i, y(x_i)) + f(x_{i+1}, y(x_{i+1})) \right] \]

带入上上式并用近似来代替式中，即可得到梯形公式

\[y_{i+1} = y_i + \frac{h}{2} \left[ f(x_i, y_i) + f(x_{i+1}, y_{i+1}) \right] \]

由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式比欧拉公式的精度高一阶

\[y_{i+1} = y_i + \frac{h}{2} \left[ f(x_i, y_i) + f(x_{i+1}, y_{i+1}) \right] \]

• 上式中的右端含有未知的\(y_{i+1}\)，它是一个关于\(y_{i+1}\)的函数方程，这类数值方法称为隐式方法。相反地，欧拉法是关于\(y_{i+1}\)的一个直接的计算公式，这类数值方法称为显示方法

两步欧拉公式

对方程\(y'=f(x,y)\)的两端在区间\([x_{i-1},x_{i+1}]\)上进行积分得到

\[y(x_{i+1}) = y(x_{i-1}) + \int_{x_{i-1}}^{x_{i+1}} f[x, y(x)] \, dx \]

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

\[\int_{x_{i-1}}^{x_{i+1}} f[x, y(x)] \, dx \approx (x_{i+1} - x_{i-1}) f[x_i, y(x_i)] \]

代入上式，并用\(y_i\)近似代替式中\(y(x_i)\)即可得到两步欧拉公式

\[y_{i+1} = y_{i-1} + 2hf(x_i, y_i) \]

前面介绍过的数值方法，无论是欧拉方法还是梯形方法，它们都是单步法，其特点是在计算\(y_{i+1}\)时只用到前一步的信息\(y_i\)，而两步欧拉公式中用到了前两步的信息

欧拉法的局部截断误差

衡量求解公式好坏的一个主要标准是求解公式的精度，因此引入局部截断误差和阶数的概念

• 对于欧拉公式，假定\(y_i=y(x_i)\)，则有

\[y_{i+1} = y(x_i) + h \left[ f(x_i, y(x_i)) \right] = y(x_i) + h y'(x_i) \]

• 而将真解\(y(x)\)在\(x_i\)处二阶泰勒展开

\[y(x_{i+1}) = y(x_i) + h y'(x_i) + \frac{h^2}{2!} y''(\xi) \quad \xi \in (x_i, x_{i+1}) \]

• 因此有

\[y(x_{i+1}) - y_{i+1} = \frac{h^2}{2!} y''(\xi) \]

• 两步欧拉公式比欧拉公式精度高一阶，设\(y_i=y(x),y_{i-1}=y(x_{i-1})\)前两步准确，则两步欧拉公式

\[y_{i+1} = y(x_{i-1}) + 2h f(x_i, y(x_i)) \]

• 把\(y(x_{i-1})\)在\(x_i\)处展开成泰勒级数，即

\[y(x_{i-1}) = y(x_i) + y'(x_i)(x_{i-1} - x_i) + \frac{y''(x_i)}{2!}(x_{i-1} - x_i)^2 + O(h^3) \]

\[\begin{align*} y(x_{i-1}) &= y(x_i) + y'(x_i)(x_{i-1} - x_i) + \frac{y''(x_i)}{2!}(x_{i-1} - x_i)^2 + O(h^3) \\ &= y(x_i) - hy'(x_i) + \frac{h^2}{2!} y''(x_i) + O(h^3) \end{align*}\]

• 由\(y_{i+1} = y(x_{i-1}) + 2hf(x_i, y_i)\)

\[\begin{align*} &= y(x_i) - hy'(x_i) + \frac{h^2}{2!}y''(x_i) + O(h^3) + 2hf(x_i, y_i) \\ &= y(x_i) + hy'(x_i) + \frac{h^2}{2!}y''(x_i) + O(h^3) \end{align*}\]

• 再将\(y(x{i+1})\)在\(x_{i+1}\)处泰勒展开

\[y(x_{i+1}) = y(x_i) + hy'(x_i) + \frac{h^2}{2!}y''(x_i) + O(h^3) \]

• 所以，可得两步欧拉公式的局部截断误差

\[y(x_{i+1}) - y_{i+1} = O(h^3) \]

改进的欧拉公式

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

• 先用欧拉公式求出一个初步的近似值\(\overline{y_{i+1}}\),称为预测值，它的精度不高，再用梯形公式对它矫正一次，即迭代一次，求得\(y_{i+1}\)称为校正值，这种预测-矫正方法称为改进的欧拉公式
  
• 可以证明，这种方法的精度为二阶。这是一种显式格式，它可以表示为嵌套形式

\[y_{i+1} = y_i + \frac{h}{2} \left[ f(x_i, y_i) + f(x_{i+1}, y_i + hf(x_i, y_i)) \right] \]

• 或者表示成下列平均化形式

\[\left\{ \begin{array}{l} y_p = y_i + h f(x_i, y_i) \\ y_c = y_i + h f(x_{i+1}, y_p) \\ y_{i+1} = \frac{1}{2}(y_p + y_c) \end{array} \right.\]

例题：用该进欧拉法解初值问题