昆仑山:眼中无形心中有穴之穴人合一

夫君子之行，静以修身，俭以养德；非澹泊无以明志，非宁静无以致远。夫学须静也，才须学也；非学无以广才，非志无以成学。怠慢则不能励精，险躁则不能冶性。年与时驰，意与岁去，遂成枯落，多不接世。悲守穷庐，将复何及！

 

9.5线性多步法

线性多步法 详细讲解与推导证明

各位同学，今天我们来系统讲解常微分方程数值解法中的线性多步法。我会从基本思想出发，完成核心公式的完整推导，明确精度、相容性的核心判定条件，最后用表格归纳所有核心知识点，帮大家把这个知识点学透、用熟。

一、线性多步法的基本思想

我们要解决的核心问题是常微分方程初值问题：

\[\begin{cases} y' = f(x,y), \quad x \in [a,b] \\ y(x_0) = y_0 \end{cases} \]

我们将求解区间\([a,b]\)做等距剖分，节点为\(x_n = x_0 + nh\)（\(h\)为固定步长，\(n=0,1,2,\dots\)），\(y_n\)是精确解\(y(x_n)\)的数值近似，\(f_n = f(x_n, y_n)\)是\(y'(x_n)\)的近似。

此前我们学的欧拉法、龙格-库塔法都属于单步法：计算\(y_{n+1}\)时，只用到前一步\(y_n\)的信息。而线性多步法的核心思想是：

计算\(y_{n+k}\)时，充分利用已经算出的\(y_n, y_{n+1}, \dots, y_{n+k-1}\)这\(k\)步的信息，通过线性组合构造递推公式，在不增加函数求值次数的前提下，获得更高的计算精度。

二、线性多步法的一般公式

1. 定义与一般形式

如果计算\(y_{n+k}\)时，除了用到\(y_{n+k-1}\)，还用到\(y_{n+i}\ (i=0,1,\dots,k-2)\)的信息，就称该方法为线性\(k\)步法，其一般公式为：

\[y_{n+k} = \sum_{i=0}^{k-1} \alpha_i y_{n+i} + h \sum_{i=0}^{k} \beta_i f_{n+i}, \quad n=0,1,2,\dots \tag{9.34} \]

符号与关键概念说明

  • \(k\)：步法的步数，\(k=1\)时退化为单步法，\(k\geq2\)为多步法；
  • \(x_{n+i} = x_n + ih\)：等距节点，\(h\)为步长；
  • \(y_{n+i}\)：精确解\(y(x_{n+i})\)的数值近似；
  • \(f_{n+i} = f(x_{n+i}, y_{n+i})\)：微分方程右端在对应节点的取值，对应\(y'\)的近似；
  • \(\alpha_i, \beta_i\)：线性多步法的待定常数，要求\(\alpha_0\)和\(\beta_0\)不全为零（否则退化为更少步数的方法）。

显式与隐式的核心区分

类型 判定条件 计算特点
显式\(k\)步法 \(\beta_k = 0\) 公式右端不含\(f_{n+k}\)，\(y_{n+k}\)可由前序步的已知值直接计算，无需迭代
隐式\(k\)步法 \(\beta_k \neq 0\) 公式右端含\(f_{n+k}=f(x_{n+k}, y_{n+k})\)，未知量\(y_{n+k}\)出现在方程两侧，需迭代求解

三、局部截断误差、精度阶与相容性（核心推导）

线性多步法的核心是通过系数\(\alpha_i, \beta_i\)的设计，控制方法的精度，我们通过泰勒展开完成完整推导。

1. 局部截断误差的定义

设\(y(x)\)是初值问题的精确解，线性多步法在\(x_{n+k}\)处的局部截断误差定义为：

\[T_{n+k} = L[y(x_n); h] = y(x_{n+k}) - \sum_{i=0}^{k-1} \alpha_i y(x_{n+i}) - h \sum_{i=0}^{k} \beta_i y'(x_{n+i}) \tag{9.35} \]

这里的“局部”，是指我们假设前序步的数值解都是精确的（\(y_{n+i}=y(x_{n+i})\)，无累积误差），此时计算\(y_{n+k}\)产生的误差，仅衡量方法本身的数学精度。

2. 精度阶的定义

若局部截断误差满足\(T_{n+k}=O(h^{p+1})\)，则称该方法是\(p\)阶精度的（简称\(p\)阶方法）；
若\(p\geq1\)，则称该方法与原微分方程是相容的（相容是方法收敛的必要条件）。

3. 泰勒展开与系数推导

我们将\(T_{n+k}\)在\(x_n\)处做泰勒展开，推导精度阶的判定条件。

首先，一元函数\(y(x)\)在\(x_n\)处的泰勒展开式为：

\[y(x_n + ih) = \sum_{q=0}^{\infty} \frac{(ih)^q}{q!} y^{(q)}(x_n) \]

其中\(y^{(0)}(x_n)=y(x_n)\)，\(y^{(q)}(x_n)\)是\(y(x)\)在\(x_n\)处的\(q\)阶导数。

对\(y'(x_n + ih)\)做泰勒展开（对上述式子求导即可）：

\[y'(x_n + ih) = \sum_{q=1}^{\infty} \frac{(ih)^{q-1}}{(q-1)!} y^{(q)}(x_n) \]

将两个展开式代入局部截断误差的定义式(9.35)，分三部分展开：

  1. \(y(x_{n+k}) = y(x_n + kh) = \sum_{q=0}^{\infty} \frac{(kh)^q}{q!} y^{(q)}(x_n)\)
  2. \(\sum_{i=0}^{k-1} \alpha_i y(x_{n+i}) = \sum_{q=0}^{\infty} \left( \frac{h^q}{q!} \sum_{i=0}^{k-1} \alpha_i i^q \right) y^{(q)}(x_n)\)（交换求和顺序）
  3. \(h \sum_{i=0}^{k} \beta_i y'(x_{n+i}) = \sum_{q=1}^{\infty} \left( \frac{h^q}{(q-1)!} \sum_{i=0}^{k} \beta_i i^{q-1} \right) y^{(q)}(x_n)\)（\(h\)与\((ih)^{q-1}\)合并为\(h^q i^{q-1}\)）

将三部分代入\(T_{n+k}\)，按\(h^q y^{(q)}(x_n)\)合并同类项，得到：

\[T_{n+k} = c_0 y(x_n) + c_1 h y'(x_n) + c_2 h^2 y''(x_n) + \dots + c_p h^p y^{(p)}(x_n) + \dots \]

其中，各阶系数\(c_q\)的表达式为：

\[\begin{cases} \displaystyle c_0 = 1 - \left( \alpha_0 + \alpha_1 + \dots + \alpha_{k-1} \right), \\ \displaystyle c_1 = k - \left[ \alpha_1 + 2\alpha_2 + \dots + (k-1)\alpha_{k-1} \right] - \left( \beta_0 + \beta_1 + \dots + \beta_k \right), \\ \displaystyle c_q = \frac{1}{q!}\left[ k^q - \left( \alpha_1 + 2^q \alpha_2 + \dots + (k-1)^q \alpha_{k-1} \right) \right] - \frac{1}{(q-1)!}\left[ \beta_1 + 2^{q-1} \beta_2 + \dots + k^{q-1} \beta_k \right], \quad q=2,3,\dots \end{cases} \tag{9.36} \]

4. 精度阶与相容性的判定条件

（1）\(p\)阶精度的充要条件

要让方法达到\(p\)阶精度，即\(T_{n+k}=O(h^{p+1})\)，需要泰勒展开中\(h^0\)到\(h^p\)的项全部为0，\(h^{p+1}\)项非零，即：

\[c_0 = c_1 = \dots = c_p = 0, \quad c_{p+1} \neq 0 \]

此时，局部截断误差可写为：

\[T_{n+k} = c_{p+1} h^{p+1} y^{(p+1)}(x_n) + O(h^{p+2}) \tag{9.37} \]

其中，第一项\(c_{p+1} h^{p+1} y^{(p+1)}(x_n)\)称为局部截断误差主项，\(c_{p+1}\)称为误差常数。

（2）相容性的充要条件

相容性要求\(p\geq1\)，即至少一阶精度，因此需要\(c_0=0\)且\(c_1=0\)，代入系数表达式得：

\[\begin{cases} \displaystyle \alpha_0 + \alpha_1 + \dots + \alpha_{k-1} = 1, \\ \displaystyle \sum_{i=1}^{k-1} i \alpha_i + \sum_{i=0}^{k} \beta_i = k. \end{cases} \tag{9.38} \]

这是线性多步法最核心的判定式：一个线性多步法与原微分方程相容，当且仅当它的系数满足(9.38)式。

四、经典案例验证（单步法特例）

当\(k=1\)时，线性多步法退化为线性单步法，我们用欧拉法和梯形法验证上述推导的正确性。

\(k=1\)时，线性1步法的一般公式为：

\[y_{n+1} = \alpha_0 y_n + h \left( \beta_0 f_n + \beta_1 f_{n+1} \right) \]

案例1：向前欧拉法（显式单步法）

令\(\beta_1=0\)（显式），用相容性条件求系数：

  1. \(c_0=0 \implies \alpha_0=1\)；
  2. \(c_1=0 \implies 1 - (\beta_0 + 0) = 0 \implies \beta_0=1\)。

得到向前欧拉法公式：

\[y_{n+1} = y_n + h f_n \]

计算精度阶：代入系数求\(c_2\)，得\(c_2=\frac{1}{2} \neq 0\)，因此\(c_0=c_1=0, c_2\neq0\)，方法为一阶精度，局部截断误差为：

\[T_{n+1} = \frac{1}{2} h^2 y''(x_n) + O(h^3) \]

与欧拉法的经典结论完全一致。

案例2：梯形法（隐式单步法）

令\(\beta_1\neq0\)（隐式），要求二阶精度，即\(c_0=c_1=c_2=0\)，求系数：

  1. \(c_0=0 \implies \alpha_0=1\)；
  2. \(c_1=0 \implies 1 - (\beta_0 + \beta_1) = 0 \implies \beta_0 + \beta_1=1\)；
  3. \(c_2=0 \implies \frac{1}{2} - \beta_1 = 0 \implies \beta_1=\frac{1}{2}\)，因此\(\beta_0=\frac{1}{2}\)。

得到梯形法公式：

\[y_{n+1} = y_n + \frac{h}{2}(f_n + f_{n+1}) \]

计算精度阶：代入系数求\(c_3\)，得\(c_3=-\frac{1}{12} \neq 0\)，因此\(c_0=c_1=c_2=0, c_3\neq0\)，方法为二阶精度，局部截断误差为：

\[T_{n+1} = -\frac{1}{12} h^3 y'''(x_n) + O(h^4) \]

与梯形法的经典结论完全一致。

五、核心知识点归纳总结

表1：线性多步法核心概念与定义

核心概念 数学定义/公式 关键说明
线性\(k\)步法一般公式 \(\displaystyle y_{n+k} = \sum_{i=0}^{k-1} \alpha_i y_{n+i} + h \sum_{i=0}^{k} \beta_i f_{n+i}\) \(x_{n+i}=x_n+ih\)，\(f_{n+i}=f(x_{n+i},y_{n+i})\)；\(\alpha_i,\beta_i\)为常数，\(\alpha_0,\beta_0\)不全为0
显式\(k\)步法 \(\beta_k=0\) 公式右端不含\(f_{n+k}\)，\(y_{n+k}\)可直接计算，无需迭代
隐式\(k\)步法 \(\beta_k \neq 0\) 公式右端含\(f_{n+k}=f(x_{n+k},y_{n+k})\)，需迭代求解\(y_{n+k}\)
局部截断误差（LTE） \(\displaystyle T_{n+k} = y(x_{n+k}) - \sum_{i=0}^{k-1} \alpha_i y(x_{n+i}) - h \sum_{i=0}^{k} \beta_i y'(x_{n+i})\) 假设前序步数值解均为精确值时，计算\(y_{n+k}\)的误差，衡量方法本身的精度
\(p\)阶精度方法 \(T_{n+k}=O(h^{p+1})\)，即\(c_0=c_1=\dots=c_p=0, c_{p+1}\neq0\) 局部截断误差主项为\(c_{p+1}h^{p+1}y^{(p+1)}(x_n)\)，\(c_{p+1}\)为误差常数
相容性 \(p\geq1\)，即\(c_0=c_1=0\) 方法与微分方程相容的充要条件：\(\sum_{i=0}^{k-1}\alpha_i=1\)，\(\sum_{i=1}^{k-1}i\alpha_i + \sum_{i=0}^k\beta_i =k\)

表2：泰勒展开系数\(c_q\)与精度条件

系数 表达式 对应精度条件
\(c_0\) \(\displaystyle 1 - \sum_{i=0}^{k-1}\alpha_i\) \(c_0=0\)：相容性的必要条件，保证方法与微分方程的一致性
\(c_1\) \(\displaystyle k - \sum_{i=1}^{k-1}i\alpha_i - \sum_{i=0}^k\beta_i\) \(c_1=0\)：相容性的必要条件，\(c_0=c_1=0\)等价于方法相容
\(c_q\ (q\geq2)\) \(\displaystyle \frac{1}{q!}\left[k^q - \sum_{i=1}^{k-1}i^q\alpha_i\right] - \frac{1}{(q-1)!}\left[\sum_{i=1}^k i^{q-1}\beta_i\right]\) \(c_0=c_1=\dots=c_p=0, c_{p+1}\neq0\) → 方法为\(p\)阶精度

表3：经典单步法参数与精度总结

方法 步数\(k\) 显/隐式 系数\(\alpha_0,\beta_0,\beta_1\) 递推公式 精度阶\(p\) 局部截断误差主项
向前欧拉法 1 显式 \(\alpha_0=1, \beta_0=1, \beta_1=0\) \(y_{n+1}=y_n + h f_n\) 1 \(\displaystyle \frac{1}{2}h^2 y''(x_n)\)
梯形法 1 隐式 \(\alpha_0=1, \beta_0=1/2, \beta_1=1/2\) \(\displaystyle y_{n+1}=y_n + \frac{h}{2}(f_n + f_{n+1})\) 2 \(\displaystyle -\frac{1}{12}h^3 y'''(x_n)\)

补充说明

  1. 线性多步法的核心优势：通过复用前序多步的计算结果，在不增加单步函数求值次数的前提下，获得比单步法更高的精度，计算效率显著提升；
  2. 多步法的使用局限：需要\(k\)个初始值\(y_0,y_1,\dots,y_{k-1}\)，其中\(y_0\)是给定初值，其余初始值需要用同阶的单步法（如龙格-库塔法）计算，才能保证整体精度；
  3. 显式与隐式的取舍：显式方法计算简单，但数值稳定性较差；隐式方法稳定性更好，但需要迭代求解，单步计算量更大。工程中常结合两者，用显式方法做预测、隐式方法做校正，构造预测-校正法。

阿当姆斯（Adams）显式与隐式公式 详细讲解与推导

各位同学，我们承接上一节线性多步法的核心理论，今天讲解工程中最常用的线性多步法——阿当姆斯（Adams）方法，包括显式（Adams-Bashforth）和隐式（Adams-Moulton）公式。我会从方法的核心思想出发，完成完整的系数推导，解读公式的精度特性，最后对比显式与隐式方法的核心差异。

一、阿当姆斯方法的核心思想与一般形式

阿当姆斯方法是线性多步法的特殊简化形式，它的构造完全源于微分方程的积分本质，我们先从源头讲起。

对于常微分方程初值问题 \(y'=f(x,y)\)，我们在区间 \([x_{n+k-1}, x_{n+k}]\) 两端对等式积分：

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

左边积分直接得到 \(y(x_{n+k}) - y(x_{n+k-1})\)，因此有：

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

阿当姆斯方法的核心，就是用多项式插值近似被积函数 \(f(x,y(x))\)，再对插值多项式积分，得到数值递推公式。对应到线性多步法的一般形式，它的系数满足：

\[\alpha_0=\alpha_1=\dots=\alpha_{k-2}=0,\quad \alpha_{k-1}=1 \]

因此阿当姆斯方法的一般公式为：

\[y_{n+k} = y_{n+k-1} + h \sum_{i=0}^{k} \beta_i f_{n+i}, \quad n=0,1,2,\dots \tag{9.39} \]

显式与隐式的核心区分

类型 判定条件 插值节点选择 计算特点
阿当姆斯显式公式（Adams-Bashforth） \(\beta_k=0\) 仅用已计算的节点 \(x_n,x_{n+1},\dots,x_{n+k-1}\) 构造插值多项式 公式右端不含未知的 \(f_{n+k}\)，\(y_{n+k}\) 可直接计算，无需迭代
阿当姆斯隐式公式（Adams-Moulton） \(\beta_k\neq0\) 用到未知节点 \(x_{n+k}\) 构造插值多项式 公式右端含未知的 \(f_{n+k}=f(x_{n+k},y_{n+k})\)，需迭代求解 \(y_{n+k}\)

二、阿当姆斯公式的完整系数推导（以k=3为例）

我们以3步阿当姆斯方法为例，结合上一节的精度阶条件，完整推导显式与隐式公式的系数。

步骤1：确定精度条件的方程组

3步阿当姆斯方法的一般形式为：

\[y_{n+3} = y_{n+2} + h\left( \beta_0 f_n + \beta_1 f_{n+1} + \beta_2 f_{n+2} + \beta_3 f_{n+3} \right) \]

根据上一节的系数公式(9.36)，由于 \(\alpha_0=\alpha_1=0,\alpha_2=1\)，\(c_0=1-(\alpha_0+\alpha_1+\alpha_2)=0\) 自动满足。要让方法达到最高精度，需令 \(c_1=c_2=c_3=c_4=0\)，代入系数公式得到方程组：

\[\begin{cases} \beta_0 + \beta_1 + \beta_2 + \beta_3 = 1, \\ 2(\beta_1 + 2\beta_2 + 3\beta_3) = 5, \\ 3(\beta_1 + 4\beta_2 + 9\beta_3) = 19, \\ 4(\beta_1 + 8\beta_2 + 27\beta_3) = 65. \end{cases} \]

步骤2：推导3步阿当姆斯显式公式

显式公式要求 \(\beta_3=0\)（\(\beta_k=0\)），因此取前3个方程求解 \(\beta_0,\beta_1,\beta_2\)：

  1. 代入 \(\beta_3=0\)，方程组简化为：

    \[\begin{cases} \beta_0 + \beta_1 + \beta_2 = 1, \\ \beta_1 + 2\beta_2 = \frac{5}{2}, \\ \beta_1 + 4\beta_2 = \frac{19}{3}. \end{cases} \]

  2. 消元求解：
    • 第3式减第2式，得 \(2\beta_2 = \frac{23}{6} \implies \beta_2 = \frac{23}{12}\)；
    • 代入第2式，得 \(\beta_1 = -\frac{16}{12}\)；
    • 代入第1式，得 \(\beta_0 = \frac{5}{12}\)。

将系数代入一般形式，得到3步阿当姆斯显式公式：

\[y_{n+3} = y_{n+2} + \frac{h}{12}\left( 23f_{n+2} - 16f_{n+1} + 5f_n \right), \quad n=0,1,2,\dots \tag{9.40} \]

精度与局部截断误差

我们令 \(c_1=c_2=c_3=0\)，代入系数计算得 \(c_4=\frac{3}{8}\neq0\)，因此该方法为三阶精度，局部截断误差为：

\[T_{n+3} = \frac{3}{8}h^4 y^{(4)}(x_n) + O(h^5) \]

步骤3：推导3步阿当姆斯隐式公式

隐式公式允许 \(\beta_3\neq0\)，因此用完整的4个方程求解4个系数：

  1. 消元求解方程组：
    • 第3式减第2式，得 \(\beta_2 + 3\beta_3 = \frac{23}{12}\)；
    • 第4式减第3式，得 \(2\beta_2 + 9\beta_3 = \frac{119}{24}\)；
    • 联立解得 \(\beta_3=\frac{3}{8}\)，\(\beta_2=\frac{19}{24}\)；
    • 回代得 \(\beta_1=-\frac{5}{24}\)，\(\beta_0=\frac{1}{24}\)。

将系数代入一般形式，得到3步阿当姆斯隐式公式：

\[y_{n+3} = y_{n+2} + \frac{h}{24}\left( 9f_{n+3} + 19f_{n+2} - 5f_{n+1} + f_n \right), \quad n=0,1,2,\dots \]

精度与局部截断误差

我们令 \(c_1=c_2=c_3=c_4=0\)，代入系数计算得 \(c_5=-\frac{19}{720}\neq0\)，因此该方法为四阶精度，局部截断误差为：

\[T_{n+3} = -\frac{19}{720}h^5 y^{(5)}(x_n) + O(h^6) \]

三、阿当姆斯公式通用表解读

我们用同样的方法，可以推导出k=1,2,3,4时的阿当姆斯显式与隐式公式，对应教材中的表9-6和表9-7，下面做完整解读。

表9-6 阿当姆斯显式公式（Adams-Bashforth）

步数k 精度阶p 递推公式 误差常数\(c_{p+1}\)
1 1 \(y_{n+1}=y_n + h f_n\)（向前欧拉法） \(\frac{1}{2}\)
2 2 \(y_{n+2}=y_{n+1} + \frac{h}{2}(3f_{n+1}-f_n)\) \(\frac{5}{12}\)
3 3 \(y_{n+3}=y_{n+2} + \frac{h}{12}(23f_{n+2}-16f_{n+1}+5f_n)\) \(\frac{3}{8}\)
4 4 \(y_{n+4}=y_{n+3} + \frac{h}{24}(55f_{n+3}-59f_{n+2}+37f_{n+1}-9f_n)\) \(\frac{251}{720}\)

显式公式核心规律：k步显式阿当姆斯方法的精度阶p=k，步数与精度阶相等。

表9-7 阿当姆斯隐式公式（Adams-Moulton）

步数k 精度阶p 递推公式 误差常数\(c_{p+1}\)
1 2 \(y_{n+1}=y_n + \frac{h}{2}(f_{n+1}+f_n)\)（梯形法） \(-\frac{1}{12}\)
2 3 \(y_{n+2}=y_{n+1} + \frac{h}{12}(5f_{n+2}+8f_{n+1}-f_n)\) \(-\frac{1}{24}\)
3 4 \(y_{n+3}=y_{n+2} + \frac{h}{24}(9f_{n+3}+19f_{n+2}-5f_{n+1}+f_n)\) \(-\frac{19}{720}\)
4 5 \(y_{n+4}=y_{n+3} + \frac{h}{720}(251f_{n+4}+646f_{n+3}-264f_{n+2}+106f_{n+1}-19f_n)\) \(-\frac{3}{160}\)

隐式公式核心规律：k步隐式阿当姆斯方法的精度阶p=k+1，同步数下，隐式比显式高一阶精度。

四、阿当姆斯显式与隐式方法核心对比

对比维度 阿当姆斯显式方法 阿当姆斯隐式方法
精度特性 同步数下精度更低，k步对应k阶 同步数下精度更高，k步对应k+1阶
计算复杂度 无需求解方程，直接递推，计算量小 需迭代求解，单步计算量更大
数值稳定性 稳定性较差，对步长h限制严格 稳定性远优于显式方法，适合刚性方程求解
局部截断误差 误差常数绝对值更大，同阶下误差更大 误差常数绝对值更小，同阶下精度更高
工程应用 常作为预测步，为隐式方法提供初值 常作为校正步，对显式结果做精度修正

补充说明

  1. 阿当姆斯方法的核心优势：结构简单，仅对\(f(x,y)\)做线性组合，单步函数求值次数极少，在长区间求解时计算效率远高于同阶龙格-库塔法；
  2. 初始值要求：k步阿当姆斯方法需要k个初始值\(y_0,y_1,\dots,y_{k-1}\)，其中\(y_0\)是给定初值，其余初始值必须用同阶或更高阶的单步法（如龙格-库塔法）计算，否则会损失整体精度；
  3. 工程常用方案：实际计算中极少单独使用显式或隐式方法，而是将两者结合，构造阿当姆斯预测-校正法：用显式公式预测\(y_{n+k}\)的初值，再用隐式公式做1-2次校正，兼顾计算效率与精度、稳定性。

例9.6 四阶阿当姆斯方法求解初值问题 完整解析

一、题目核心信息

1. 初值问题

求解常微分方程初值问题：

\[\begin{cases} y' = -y + x + 1 \\ y(0) = 1 \end{cases} \]

取步长 \(h=0.1\)，该方程的精确解为 \(y(x) = e^{-x} + x\)。

2. 核心准备

  • 微分方程右端项：\(f_n = f(x_n,y_n) = -y_n + x_n + 1\)
  • 等距节点：\(x_n = nh = 0.1n\)
  • 求解方法：四阶阿当姆斯显式公式（Adams-Bashforth）、四阶阿当姆斯隐式公式（Adams-Moulton）

二、公式化简与推导

1. 四阶阿当姆斯显式公式化简

四阶阿当姆斯显式公式的通用形式为：

\[y_{n+4} = y_{n+3} + \frac{h}{24}\left(55f_{n+3} - 59f_{n+2} + 37f_{n+1} - 9f_n\right) \]

将 \(f_{n+i} = -y_{n+i} + x_{n+i} + 1\)、\(h=0.1\)、\(x_{n+i}=0.1(n+i)\) 代入，展开并合并同类项：

  1. 展开所有 \(f\) 项，合并 \(y\) 的线性项、\(x\) 的线性项和常数项；
  2. 代入 \(h=0.1\) 化简系数，最终得到显式递推公式：

\[y_{n+4} = \frac{1}{24}\left(18.5y_{n+3} + 5.9y_{n+2} - 3.7y_{n+1} + 0.9y_n + 0.24n + 3.24\right) \]

该公式为显式递推：已知前4个值 \(y_n,y_{n+1},y_{n+2},y_{n+3}\)，可直接计算 \(y_{n+4}\)，无需迭代。

2. 四阶阿当姆斯隐式公式化简

四阶阿当姆斯隐式公式的通用形式为：

\[y_{n+3} = y_{n+2} + \frac{h}{24}\left(9f_{n+3} + 19f_{n+2} - 5f_{n+1} + f_n\right) \]

同样代入 \(f_{n+i}\)、\(h=0.1\)、\(x_{n+i}=0.1(n+i)\) 展开：

  1. 展开后，未知量 \(y_{n+3}\) 同时出现在公式左右两侧；
  2. 由于原微分方程是线性微分方程，\(f\) 对 \(y\) 是线性的，因此可直接移项合并，解出 \(y_{n+3}\)，无需迭代，最终得到：

\[y_{n+3} = \frac{1}{24.9}\left(22.1y_{n+2} + 0.5y_{n+1} - 0.1y_n + 0.24n + 3\right) \]

若原方程为非线性微分方程，隐式公式无法直接求解，需用迭代法计算 \(y_{n+3}\)。

三、初始值计算

线性多步法需要前置初始值才能启动递推：

  • 四阶显式方法为4步法，需要4个初始值 \(y_0,y_1,y_2,y_3\)；
  • 四阶隐式方法为3步法，需要3个初始值 \(y_0,y_1,y_2\)。

本题有精确解，直接代入精确解 \(y(x)=e^{-x}+x\) 计算初始值：

\(x_n\) 0.0 0.1 0.2 0.3
\(y_n\) 1.0 \(e^{-0.1}+0.1\approx1.004837418\) \(e^{-0.2}+0.2\approx1.018730753\) \(e^{-0.3}+0.3\approx1.040818221\)

实际工程中若无精确解，需用同阶单步法（如四阶龙格-库塔法）计算初始值，保证初始值精度与多步法匹配。

四、计算结果与误差分析

1. 完整计算结果表

\(x_n\) 精确解 \(y(x_n)=e^{-x_n}+x_n\) 阿当姆斯显式方法 阿当姆斯隐式方法
\(y_n\) 绝对误差 \(|y(x_n)-y_n|\) \(y_n\) 绝对误差 \(|y(x_n)-y_n|\)
0.3 1.040 818 22 - - 1.040 818 01 \(2.1\times10^{-7}\)
0.4 1.070 320 05 1.070 322 92 \(2.87\times10^{-6}\) 1.070 319 66 \(3.8\times10^{-7}\)
0.5 1.106 530 66 1.106 535 48 \(4.82\times10^{-6}\) 1.106 530 14 \(5.2\times10^{-7}\)
0.6 1.148 811 64 1.148 818 41 \(6.77\times10^{-6}\) 1.148 811 01 \(6.3\times10^{-7}\)
0.7 1.196 585 30 1.196 593 39 \(8.09\times10^{-6}\) 1.196 584 59 \(7.1\times10^{-7}\)
0.8 1.249 328 96 1.249 338 16 \(9.19\times10^{-6}\) 1.249 328 19 \(7.7\times10^{-7}\)
0.9 1.306 569 66 1.306 579 61 \(9.95\times10^{-6}\) 1.306 568 85 \(8.1\times10^{-7}\)
1.0 1.367 879 44 1.367 889 96 \(1.05\times10^{-5}\) 1.367 878 60 \(8.4\times10^{-7}\)

2. 核心结果分析

  1. 精度差异：同是四阶精度，阿当姆斯隐式方法的绝对误差比显式方法小一个数量级（隐式误差在 \(10^{-7}\) 量级，显式在 \(10^{-6}\sim10^{-5}\) 量级）。
  2. 误差来源：两者局部截断误差主项的系数差异显著：
    • 四阶显式阿当姆斯的误差常数 \(c_{p+1}=\frac{251}{720}\approx0.3486\)
    • 四阶隐式阿当姆斯的误差常数 \(c_{p+1}=-\frac{19}{720}\approx-0.0264\)
      隐式方法的误差常数绝对值远小于显式，因此局部截断误差更小，整体精度更高。
  3. 计算效率：本题为线性微分方程，隐式公式可直接求解，无需迭代，计算量与显式相当，但精度优势显著；若为非线性方程，隐式方法需迭代求解，单步计算量会大于显式方法。

五、核心知识点总结

  1. 多步法的启动条件：线性多步法无法自启动，必须用同阶单步法（如龙格-库塔法）计算前置初始值，初始值的精度直接影响整体计算结果。
  2. 显式与隐式的取舍：显式方法计算简单、无需迭代，但稳定性差、精度偏低；隐式方法精度高、稳定性好，线性方程可直接求解，非线性方程需迭代。
  3. 工程常用方案：实际计算中常采用预测-校正法：用显式公式预测 \(y_{n+k}\) 的初值，再用隐式公式做1~2次校正，兼顾计算效率与精度、稳定性。

米尔尼方法与辛普森方法 详细讲解与推导

本节我们讲解两种经典的线性多步法：米尔尼（Milne）方法和辛普森（Simpson）方法，二者均通过“微分方程积分+数值求积/泰勒展开”构造，是常微分方程初值问题数值求解的常用方法，我们将从构造思路、完整推导、精度特性、核心特点四个维度展开讲解。

一、米尔尼（Milne）方法

米尔尼方法是四阶四步显式线性多步法，是阿当姆斯显式方法之外最常用的显式多步法，我们通过泰勒展开系数匹配法完成完整推导。

1. 方法的构造形式

米尔尼方法是k=4的显式线性多步法，其构造形式为：

\[y_{n+4} = y_n + h\left( \beta_3 f_{n+3} + \beta_2 f_{n+2} + \beta_1 f_{n+1} + \beta_0 f_n \right) \]

对应线性多步法的一般形式，其系数满足：\(\alpha_0=1,\ \alpha_1=\alpha_2=\alpha_3=0\)，\(\beta_4=0\)（显式，不含\(f_{n+4}\)）。我们的目标是确定\(\beta_0,\beta_1,\beta_2,\beta_3\)，使方法的精度阶尽可能高。

2. 系数的完整推导

回顾线性多步法的精度阶判定条件：要使方法达到p阶精度，需满足局部截断误差的泰勒展开系数\(c_0=c_1=\dots=c_p=0\)，\(c_{p+1}\neq0\)。

首先代入\(\alpha_0=1,\alpha_1=\alpha_2=\alpha_3=0\)，计算各阶系数：

  1. \(c_0\)的计算：

    \[c_0 = 1 - \sum_{i=0}^{3}\alpha_i = 1 - 1 = 0 \]

    自动满足相容性的必要条件。

  2. 令\(c_1=c_2=c_3=c_4=0\)，建立方程组
    结合线性多步法的系数公式(9.36)，代入k=4和\(\alpha_i\)的取值，得到：

    \[\begin{cases} \displaystyle c_1 = 4 - \sum_{i=0}^{3}\beta_i = 0 \implies \beta_0 + \beta_1 + \beta_2 + \beta_3 = 4, \\ \displaystyle c_2 = \frac{4^2}{2!} - \sum_{i=1}^{3}i\beta_i = 0 \implies 2\left( \beta_1 + 2\beta_2 + 3\beta_3 \right) = 16, \\ \displaystyle c_3 = \frac{4^3}{3!} - \frac{1}{2!}\sum_{i=1}^{3}i^2\beta_i = 0 \implies 3\left( \beta_1 + 4\beta_2 + 9\beta_3 \right) = 64, \\ \displaystyle c_4 = \frac{4^4}{4!} - \frac{1}{3!}\sum_{i=1}^{3}i^3\beta_i = 0 \implies 4\left( \beta_1 + 8\beta_2 + 27\beta_3 \right) = 256. \end{cases} \]

  3. 解线性方程组
    采用消元法求解：

    • 化简后三个方程，消去\(\beta_1\)，得到关于\(\beta_2,\beta_3\)的方程组，解得\(\beta_3=\frac{8}{3}\)，\(\beta_2=-\frac{4}{3}\)；
    • 回代得\(\beta_1=\frac{8}{3}\)，再代入第一个方程得\(\beta_0=0\)。

    最终系数为：

    \[\beta_0=0,\quad \beta_1=\frac{8}{3},\quad \beta_2=-\frac{4}{3},\quad \beta_3=\frac{8}{3} \]

3. 米尔尼方法的最终公式与精度

将系数代入构造形式，得到米尔尼方法的递推公式：

\[y_{n+4} = y_n + \frac{4h}{3}\left( 2f_{n+3} - f_{n+2} + 2f_{n+1} \right), \quad n=0,1,2,\dots \tag{9.41} \]

精度与局部截断误差

计算\(c_5\)（\(c_0\)到\(c_4\)均为0）：

\[c_5 = \frac{4^5}{5!} - \frac{1}{4!}\sum_{i=1}^{3}i^4\beta_i = \frac{1024}{120} - \frac{1}{24}\left( \frac{8}{3} - 16\times\frac{4}{3} + 81\times\frac{8}{3} \right) = \frac{14}{45} \neq 0 \]

因此米尔尼方法为四阶精度，局部截断误差为：

\[T_{n+4} = \frac{14}{45}h^5 y^{(5)}(x_n) + O(h^6) \tag{9.42} \]

4. 积分构造思路

米尔尼方法也可通过微分方程的积分直接构造：对\(y'=f(x,y)\)在区间\([x_n, x_{n+4}]\)两端积分，得

\[y(x_{n+4}) - y(x_n) = \int_{x_n}^{x_{n+4}} f(x,y(x))dx \]

对右端被积函数用三次插值多项式近似，积分后即可得到与泰勒展开法完全一致的米尔尼公式。

二、辛普森（Simpson）方法

辛普森方法是四阶二步隐式线性多步法，直接源于数值积分中的辛普森求积公式，是精度极高的隐式多步法。

1. 方法的构造思路

对常微分方程\(y'=f(x,y)\)在区间\([x_n, x_{n+2}]\)两端积分，得：

\[y(x_{n+2}) - y(x_n) = \int_{x_n}^{x_{n+2}} f(x,y(x))dx \]

对右端的定积分，采用辛普森求积公式近似。辛普森求积公式是二阶牛顿-科茨公式，形式为：

\[\int_{a}^{b} g(x)dx \approx \frac{b-a}{6}\left[ g(a) + 4g\left( \frac{a+b}{2} \right) + g(b) \right] \]

此处\(a=x_n\)，\(b=x_{n+2}\)，区间长度\(b-a=2h\)，中点为\(x_{n+1}\)，代入得：

\[\int_{x_n}^{x_{n+2}} f(x,y(x))dx \approx \frac{2h}{6}\left[ f(x_n,y(x_n)) + 4f(x_{n+1},y(x_{n+1})) + f(x_{n+2},y(x_{n+2})) \right] \]

2. 辛普森方法的最终公式与精度

将积分近似代入原式，得到辛普森方法的递推公式：

\[y_{n+2} = y_n + \frac{h}{3}\left( f_n + 4f_{n+1} + f_{n+2} \right), \quad n=0,1,2,\dots \]

方法类型与精度

  • 方法类型：k=2步隐式线性多步法（公式右端含未知的\(f_{n+2}=f(x_{n+2},y_{n+2})\)）；
  • 精度特性：辛普森求积公式的截断误差为\(-\frac{(b-a)}{180}\left( \frac{b-a}{2} \right)^4 g^{(4)}(\xi)\)，代入\(b-a=2h\)，可得局部截断误差为：

    \[T_{n+2} = -\frac{h^5}{90}y^{(5)}(x_n) + O(h^6) \]

    因此辛普森方法为四阶精度，是二步方法能达到的最高精度。

三、两种方法的核心特性对比与总结

表1 米尔尼方法与辛普森方法核心参数对比

特性 米尔尼方法 辛普森方法
步数k 4步 2步
显/隐式 显式（\(\beta_4=0\)） 隐式（\(\beta_2\neq0\)）
精度阶p 四阶 四阶
递推公式 \(\displaystyle y_{n+4}=y_n+\frac{4h}{3}(2f_{n+3}-f_{n+2}+2f_{n+1})\) \(\displaystyle y_{n+2}=y_n+\frac{h}{3}(f_n+4f_{n+1}+f_{n+2})\)
局部截断误差主项 \(\displaystyle \frac{14}{45}h^5 y^{(5)}(x_n)\) \(\displaystyle -\frac{h^5}{90}y^{(5)}(x_n)\)
误差常数绝对值 \(\approx0.311\) \(\approx0.011\)
计算特点 直接递推，无需迭代，计算量小 需迭代求解，线性方程可直接化简求解
初始值需求 需4个初始值\(y_0,y_1,y_2,y_3\) 需2个初始值\(y_0,y_1\)

核心结论

  1. 精度差异：二者同为四阶精度，但辛普森方法的误差常数绝对值远小于米尔尼方法，同步长下辛普森方法的局部截断误差更小，精度更高；
  2. 计算效率：米尔尼方法为显式，无需迭代，单步计算量更小，但需要更多初始值；辛普森方法为隐式，非线性方程需迭代，但初始值需求少，数值稳定性更优；
  3. 工程应用：米尔尼方法常作为预测步，与汉明（Hamming）校正法结合，构成米尔尼-汉明预测-校正法；辛普森方法是隐式多步法的基础，也是刚性微分方程求解的常用方法，同时是数值积分与微分方程数值解的核心桥梁。

汉明（Hamming）方法 详细讲解与推导

汉明方法是为解决辛普森方法精度高但数值稳定性差的问题而提出的经典线性多步法，它是三步四阶隐式线性多步法，在保持四阶高精度的同时大幅提升了数值稳定性，是常微分方程初值问题数值求解中最常用的多步法之一。

一、方法的构造背景与一般形式

1. 提出背景

辛普森方法是二步四阶方法，是同步数下精度最高的方法，但存在数值稳定性差的缺陷，在长区间、大步长求解时容易出现数值发散。为改善稳定性，我们构造一类三步线性多步法，在保证四阶精度的前提下优化稳定性。

2. 一般构造形式

三步线性多步法的一般形式为：

\[y_{n+3} = \alpha_0 y_n + \alpha_1 y_{n+1} + \alpha_2 y_{n+2} + h\left( \beta_1 f_{n+1} + \beta_2 f_{n+2} + \beta_3 f_{n+3} \right) \]

其中：

  • \(k=3\)为步数，\(f_{n+i}=f(x_{n+i},y_{n+i})\)，\(x_{n+i}=x_n+ih\)；
  • \(\alpha_0,\alpha_1,\alpha_2,\beta_1,\beta_2,\beta_3\)为待定常数；
  • 若取\(\alpha_1=1\)，可退化为辛普森方法；我们取\(\alpha_1=0\)，构造全新的三步法以优化稳定性。

二、系数的完整推导（四阶精度条件）

要使方法达到四阶精度，需满足局部截断误差的泰勒展开系数\(c_0=c_1=c_2=c_3=c_4=0\)（\(c_5\neq0\)）。结合线性多步法的通用系数公式(9.36)，代入\(k=3\)、\(\alpha_1=0\)，建立方程组并求解。

步骤1：建立精度条件方程组

根据线性多步法的系数公式，四阶精度要求\(c_0\)到\(c_4\)全为0，推导得方程组：

\[\begin{cases} \alpha_0 + \alpha_2 = 1, \\ 2\alpha_2 + \beta_1 + \beta_2 + \beta_3 = 3, \\ 4\alpha_2 + 2\left( \beta_1 + 2\beta_2 + 3\beta_3 \right) = 9, \\ 8\alpha_2 + 3\left( \beta_1 + 4\beta_2 + 9\beta_3 \right) = 27, \\ 16\alpha_2 + 4\left( \beta_1 + 8\beta_2 + 27\beta_3 \right) = 81. \end{cases} \]

步骤2：消元求解系数

  1. 由第一个方程得\(\alpha_0=1-\alpha_2\)，将剩余方程转化为关于\(\alpha_2\)和\(\beta_1,\beta_2,\beta_3\)的线性组合；
  2. 对后四个方程做消元处理，先解出\(\beta_1,\beta_2,\beta_3\)关于\(\alpha_2\)的表达式，再代入第五个方程求解\(\alpha_2\)；
  3. 最终解得所有系数：

    \[\alpha_0 = -\frac{1}{8},\quad \alpha_2 = \frac{9}{8},\quad \beta_1 = -\frac{3}{8},\quad \beta_2 = \frac{6}{8},\quad \beta_3 = \frac{3}{8} \]

三、汉明方法的最终公式与精度特性

1. 汉明方法的递推公式

将系数代入一般形式，整理得到汉明方法的最终递推公式：

\[y_{n+3} = \frac{1}{8}\left( 9y_{n+2} - y_n \right) + \frac{3h}{8}\left( f_{n+3} + 2f_{n+2} - f_{n+1} \right), \quad n=0,1,2,\dots \tag{9.43} \]

方法类型：公式右端含未知项\(f_{n+3}=f(x_{n+3},y_{n+3})\)，因此是隐式线性多步法。对于线性微分方程，可直接移项解出\(y_{n+3}\)，无需迭代；对于非线性微分方程，需迭代求解。

2. 精度与局部截断误差

计算\(c_5\)（\(c_0\)到\(c_4\)均为0），得\(c_5=-\frac{1}{40}\neq0\)，因此汉明方法为四阶精度，其局部截断误差为：

\[T_{n+3} = -\frac{h^5}{40}y^{(5)}(x_n) + O(h^6) \tag{9.44} \]

其中，\(-\frac{h^5}{40}y^{(5)}(x_n)\)为局部截断误差主项，\(-\frac{1}{40}\)为误差常数。

四、汉明方法的核心特性与工程应用

1. 核心优势

  1. 精度优秀：同为四阶精度，汉明方法的误差常数绝对值（\(1/40=0.025\)）远小于米尔尼方法（\(14/45\approx0.311\)），略优于四阶阿当姆斯隐式方法（\(19/720\approx0.026\)），同步长下精度更高；
  2. 稳定性优异：核心设计目标是改善辛普森方法的稳定性缺陷，汉明方法的数值稳定区间远大于同阶的辛普森方法、米尔尼方法，适合长区间、大步长的微分方程求解；
  3. 计算效率均衡：三步法仅需3个初始值即可启动，初始值需求少于四步的米尔尼方法、阿当姆斯方法，工程实现更便捷。

2. 典型工程应用

汉明方法最经典的应用是米尔尼-汉明预测-校正法：

  • 预测步：用四阶米尔尼显式公式计算\(y_{n+4}\)的预测初值；
  • 校正步：用汉明隐式公式对预测值做1~2次校正，得到高精度结果。

该方案结合了显式方法的计算效率和隐式方法的精度、稳定性，是工业界求解非刚性常微分方程的主流方案之一。

五、核心参数总结表

特性 汉明方法
步数\(k\) 3步
显/隐式 隐式（含\(f_{n+3}\)）
精度阶\(p\) 四阶
递推公式 \(\displaystyle y_{n+3} = \frac{1}{8}\left( 9y_{n+2} - y_n \right) + \frac{3h}{8}\left( f_{n+3} + 2f_{n+2} - f_{n+1} \right)\)
局部截断误差主项 \(\displaystyle -\frac{h^5}{40}y^{(5)}(x_n)\)
误差常数绝对值 \(0.025\)
初始值需求 3个初始值\(y_0,y_1,y_2\)
核心优势 四阶高精度、数值稳定性优异、初始值需求少

预测-校正方法 详细讲解与推导

预测-校正方法是线性多步法的工程化实现方案，核心是结合显式多步法的计算便捷性与隐式多步法的高精度、高稳定性，在避免迭代计算的同时，获得接近隐式方法的求解效果，是工业界求解常微分方程初值问题的主流方案之一。

一、方法的核心思想与基本逻辑

1. 提出背景

  • 隐式线性多步法（如阿当姆斯隐式、汉明方法）：精度高、数值稳定性好，但求解时需迭代，非线性微分方程的迭代计算量极大；
  • 显式线性多步法（如阿当姆斯显式、米尔尼方法）：无需迭代，单步计算量小，但精度、数值稳定性弱于同阶隐式方法。

预测-校正方法的核心是扬长避短：用同阶精度的显式公式给出未知量的初始近似（预测），再用隐式公式对近似值做一次修正（校正），无需迭代即可获得接近隐式方法的精度与稳定性，同时保持显式方法的计算效率。

2. 核心步骤定义

预测-校正方法的核心环节可缩写为4类，分别是：

缩写 环节名称 核心作用
P 预测Predictor 用显式多步法公式，计算\(y_{n+k}\)的初始预测值\(y_{n+k}^p\)，无需迭代
E 求值Evaluation 用预测/校正后的\(y\)值，计算微分方程右端项\(f=f(x,y)\)，为后续步骤提供输入
C 校正Corrector 用隐式多步法公式，代入预测得到的\(f\)值，计算校正后的\(y_{n+k}^c\)，仅做一次修正，无需迭代
M 修正Modify 利用预测值与校正值的差值，补偿截断误差，进一步提升精度、抑制误差累积

3. 基础原型：二阶改进欧拉法

最简单的预测-校正方法是二阶改进欧拉法，完美体现了核心逻辑：

  • 预测P（向前欧拉法，显式一阶）：\(y_{n+1}^p = y_n + h f_n\)
  • 求值E：\(f_{n+1}^p = f(x_{n+1}, y_{n+1}^p)\)
  • 校正C（梯形法，隐式二阶）：\(y_{n+1} = y_n + \frac{h}{2}(f_n + f_{n+1}^p)\)
  • 求值E：\(f_{n+1} = f(x_{n+1}, y_{n+1})\)

核心原则：预测公式与校正公式必须为同阶精度的显式-隐式匹配，否则会损失整体求解精度。

二、四阶阿当姆斯预测-校正格式（PECE）

四阶阿当姆斯PECE格式是最经典的四阶预测-校正方案，采用四阶阿当姆斯显式公式（4步4阶）做预测，四阶阿当姆斯隐式公式（3步4阶）做校正，二者精度匹配，结构简单易实现。

1. 完整计算步骤

\[\begin{align*} &\text{预测P：} \quad y_{n+4}^p = y_{n+3} + \frac{h}{24}\left(55f_{n+3} - 59f_{n+2} + 37f_{n+1} - 9f_n\right), \\ &\text{求值E：} \quad f_{n+4}^p = f\left(x_{n+4}, y_{n+4}^p\right), \\ &\text{校正C：} \quad y_{n+4} = y_{n+3} + \frac{h}{24}\left(9f_{n+4}^p + 19f_{n+3} - 5f_{n+2} + f_{n+1}\right), \\ &\text{求值E：} \quad f_{n+4} = f\left(x_{n+4}, y_{n+4}\right). \end{align*} \]

2. 事后误差估计推导

预测-校正方法的核心优势之一，是可通过预测值与校正值的差值，无需精确解即可估算局部截断误差，为步长自适应提供依据。

首先写出两步的局部截断误差主项：

  • 预测步（阿当姆斯显式）：\(\displaystyle y(x_{n+4}) - y_{n+4}^p \approx \frac{251}{720}h^5 y^{(5)}(x_n)\)
  • 校正步（阿当姆斯隐式）：\(\displaystyle y(x_{n+4}) - y_{n+4} \approx -\frac{19}{720}h^5 y^{(5)}(x_n)\)

两式相减，消去未知的精确解\(y(x_{n+4})\)，得：

\[y_{n+4} - y_{n+4}^p \approx \frac{270}{720}h^5 y^{(5)}(x_n) \]

整理得到截断误差核心项的近似：

\[h^5 y^{(5)}(x_n) \approx -\frac{720}{270}\left(y_{n+4}^p - y_{n+4}\right) \]

代回误差公式，得到事后误差估计式：

  • 预测值误差：\(\displaystyle y(x_{n+4}) - y_{n+4}^p \approx -\frac{251}{270}\left(y_{n+4}^p - y_{n+4}\right)\)
  • 校正值误差：\(\displaystyle y(x_{n+4}) - y_{n+4} \approx \frac{19}{270}\left(y_{n+4}^p - y_{n+4}\right)\)

3. 格式特点

  • 优势：单步计算量小，仅需2次函数求值，可实时估算计算误差，易实现；
  • 局限：无误差修正，误差累积较快，数值稳定性一般，长区间求解易发散。

三、修正型阿当姆斯预测-校正格式（PMECME）

基于PECE格式的误差估计，我们可以通过误差补偿修正进一步提升精度、抑制误差累积，构造带修正环节的PMECME格式。

1. 修正的核心逻辑

通过事后误差估计，我们可以直接补偿预测值和校正值的截断误差：

  • 预测值修正：\(\displaystyle y_{n+4}^{pm} = y_{n+4}^p + \frac{251}{270}\left(y_{n+4} - y_{n+4}^p\right)\)，补偿预测值的误差；
  • 校正值修正：\(\displaystyle y_{n+4} = y_{n+4}^c - \frac{19}{270}\left(y_{n+4}^c - y_{n+4}^p\right)\)，补偿校正值的误差。

由于修正预测值时，当前步的\(y_{n+4}\)未知，因此用上一步的预测-校正差值\(y_{n+3}^c - y_{n+3}^p\)代替，实现实时修正。

2. 完整计算步骤

\[\begin{align*} &\text{预测P：} \quad y_{n+4}^p = y_{n+3} + \frac{h}{24}\left(55f_{n+3} - 59f_{n+2} + 37f_{n+1} - 9f_n\right), \\ &\text{修正M：} \quad y_{n+4}^{pm} = y_{n+4}^p + \frac{251}{270}\left(y_{n+3}^c - y_{n+3}^p\right), \\ &\text{求值E：} \quad f_{n+4}^{pm} = f\left(x_{n+4}, y_{n+4}^{pm}\right), \\ &\text{校正C：} \quad y_{n+4}^c = y_{n+3} + \frac{h}{24}\left(9f_{n+4}^{pm} + 19f_{n+3} - 5f_{n+2} + f_{n+1}\right), \\ &\text{修正M：} \quad y_{n+4} = y_{n+4}^c - \frac{19}{270}\left(y_{n+4}^c - y_{n+4}^p\right), \\ &\text{求值E：} \quad f_{n+4} = f\left(x_{n+4}, y_{n+4}\right). \end{align*} \]

3. 格式优势

  • 带误差补偿，精度显著高于PECE格式；
  • 有效抑制误差累积，数值稳定性更优；
  • 仍仅需2次函数求值，计算量增幅极小。

四、米尔尼-汉明预测-校正格式

米尔尼-汉明格式是工程中最常用的四阶预测-校正方案，采用米尔尼显式公式（4步4阶）做预测，汉明隐式公式（3步4阶）做校正，核心优势是汉明方法的数值稳定性远优于阿当姆斯隐式方法，长区间求解表现更优异。

1. 误差修正的来源

两步的局部截断误差主项为：

  • 预测步（米尔尼公式）：\(\displaystyle y(x_{n+4}) - y_{n+4}^p \approx \frac{14}{45}h^5 y^{(5)}(x_n)\)
  • 校正步（汉明公式）：\(\displaystyle y(x_{n+4}) - y_{n+4}^c \approx -\frac{1}{40}h^5 y^{(5)}(x_n)\)

同前述方法推导，得到修正系数分别为\(\frac{112}{121}\)（预测修正）和\(-\frac{9}{121}\)（校正修正）。

2. 完整计算步骤

\[\begin{align*} &\text{预测P：} \quad y_{n+4}^p = y_n + \frac{4h}{3}\left(2f_{n+3} - f_{n+2} + 2f_{n+1}\right), \\ &\text{修正M：} \quad y_{n+4}^{pm} = y_{n+4}^p + \frac{112}{121}\left(y_{n+3}^c - y_{n+3}^p\right), \\ &\text{求值E：} \quad f_{n+4}^{pm} = f\left(x_{n+4}, y_{n+4}^{pm}\right), \\ &\text{校正C：} \quad y_{n+4}^c = \frac{1}{8}\left(9y_{n+3} - y_{n+1}\right) + \frac{3h}{8}\left(f_{n+4}^{pm} + 2f_{n+3} - f_{n+2}\right), \\ &\text{修正M：} \quad y_{n+4} = y_{n+4}^c - \frac{9}{121}\left(y_{n+4}^c - y_{n+4}^p\right), \\ &\text{求值E：} \quad f_{n+4} = f\left(x_{n+4}, y_{n+4}\right). \end{align*} \]

3. 核心优势

  • 汉明方法的稳定区间远大于阿当姆斯方法，长时程、大步长求解时不易发散；
  • 误差常数更小，同步长下精度更高；
  • 带误差修正，误差累积抑制效果好，是工业级仿真的首选方案之一。

五、核心格式对比与使用规范

1. 不同预测-校正格式核心特性对比

格式名称 预测公式 校正公式 精度阶 核心步骤 核心优势 局限性
二阶改进欧拉法 向前欧拉法 梯形法 二阶 PECE 结构极简，易实现 精度低，仅适合教学与简单场景
四阶阿当姆斯PECE 四阶阿当姆斯显式 四阶阿当姆斯隐式 四阶 PECE 计算量极小，可误差估计 稳定性一般，误差累积快
四阶阿当姆斯PMECME 四阶阿当姆斯显式 四阶阿当姆斯隐式 四阶 PMECME 精度高，误差累积小，计算量低 步骤稍多，稳定性仍有限
四阶米尔尼-汉明PMECME 米尔尼显式公式 汉明隐式公式 四阶 PMECME 稳定性优异，长区间表现好，精度高 需4个初始值，无法自启动

2. 适用场景

  1. 非刚性常微分方程（组）的长时程积分，如天体轨道计算、机械动力学仿真、生态系统演化模拟；
  2. 算力有限、需兼顾精度与实时性的场景，如嵌入式系统、硬件在环仿真、实时控制系统；
  3. 大规模微分方程组求解，多步法单步计算量远小于同阶龙格-库塔法，预测-校正进一步降低计算成本；
  4. 需要自适应步长的求解场景，事后误差估计可直接作为步长调整的依据。

3. 使用注意事项

  1. 精度匹配：预测公式与校正公式必须为同阶精度，否则会损失整体求解精度；
  2. 初始值要求：多步法无法自启动，必须用同阶或更高阶的单步法（如四阶龙格-库塔法）计算初始值，初始值精度直接决定整体求解精度；
  3. 步长控制：步长需控制在方法的稳定区间内，过大的步长会导致数值振荡、发散；
  4. 刚性方程不适用：预测-校正方法的稳定区间有限，无法应对刚性微分方程，此类场景优先使用向后差分法（BDF）等A-稳定方法。

多步法公式构造方法与例题详解

本节核心讲解任意形式线性多步法的通用构造方法——泰勒展开法，并通过两个典型例题完整演示推导过程，明确精度阶的判定与局部截断误差的计算方法。

一、构造多步法的两种途径核心对比

构造线性多步法有两种核心途径，二者的适用范围与灵活性差异显著：

构造途径 核心原理 适用范围 局限性
数值积分法 对微分方程\(y'=f(x,y)\)两端积分，用数值求积公式近似积分项 仅适用于可转化为积分方程的多步法（如阿当姆斯类方法） 通用性差，仅能构造特定形式的多步法
泰勒展开法 基于局部截断误差的定义，将所有\(y(x_{n+i})\)、\(y'(x_{n+i})\)在\(x_n\)处做泰勒展开，通过匹配低阶项系数为0确定待定参数 适用于任意形式的线性多步法，无论是否符合标准线性k步法形式 无本质局限性，是构造多步法的通用方法，不易出错

教材注记的核心结论：优先使用泰勒展开法构造多步法，无需套用固定系数公式，可灵活处理任意形式的多步法，是工程与考试中最常用的方法。

二、例9.7 显式二步法的参数确定与精度分析

1. 题目重述

求解初值问题\(y'=f(x,y),\ y(x_0)=y_0\)，采用显式二步法：

\[y_{n+1} = \alpha_0 y_n + \alpha_1 y_{n-1} + h\left( \beta_0 f_n + \beta_1 f_{n-1} \right) \]

其中\(f_n=f(x_n,y_n)\)，\(f_{n-1}=f(x_{n-1},y_{n-1})\)。试确定参数\(\alpha_0,\alpha_1,\beta_0,\beta_1\)，使方法的精度阶尽可能高，并求局部截断误差。

2. 核心思路：局部截断误差的泰勒展开

步骤1：写出局部截断误差的定义

局部截断误差的核心假设：前序步的数值解均为精确值，即\(y_n=y(x_n)\)，\(y_{n-1}=y(x_{n-1})\)，\(f_n=y'(x_n)\)，\(f_{n-1}=y'(x_{n-1})\)。因此方法在\(x_{n+1}\)处的局部截断误差为：

\[T_{n+1} = y(x_{n+1}) - \left[ \alpha_0 y(x_n) + \alpha_1 y(x_{n-1}) + h\left( \beta_0 y'(x_n) + \beta_1 y'(x_{n-1}) \right) \right] \]

其中\(x_{n+1}=x_n+h\)，\(x_{n-1}=x_n-h\)。

步骤2：泰勒展开所有项

将所有函数在\(x_n\)处做泰勒展开，为了匹配尽可能高的精度，展开到\(h^4\)项：

  1. \(y(x_n+h) = y(x_n) + h y'(x_n) + \frac{h^2}{2!}y''(x_n) + \frac{h^3}{3!}y'''(x_n) + \frac{h^4}{4!}y^{(4)}(x_n) + O(h^5)\)
  2. \(y(x_n-h) = y(x_n) - h y'(x_n) + \frac{h^2}{2!}y''(x_n) - \frac{h^3}{3!}y'''(x_n) + \frac{h^4}{4!}y^{(4)}(x_n) + O(h^5)\)
  3. \(y'(x_n-h) = y'(x_n) - h y''(x_n) + \frac{h^2}{2!}y'''(x_n) - \frac{h^3}{3!}y^{(4)}(x_n) + O(h^4)\)

步骤3：代入并合并同类项

将展开式代入\(T_{n+1}\)，按\(h\)的幂次整理，合并同类项：

\[\begin{align*} T_{n+1} &= \left[ y(x_n) + h y' + \frac{h^2}{2}y'' + \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} \right] \\ &- \alpha_0 y(x_n) - \alpha_1 \left[ y(x_n) - h y' + \frac{h^2}{2}y'' - \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} \right] \\ &- h\beta_0 y' - h\beta_1 \left[ y' - h y'' + \frac{h^2}{2}y''' - \frac{h^3}{6}y^{(4)} \right] + O(h^5) \\ &= \underbrace{(1 - \alpha_0 - \alpha_1)}_{c_0} y(x_n) \\ &+ \underbrace{(1 + \alpha_1 - \beta_0 - \beta_1)}_{c_1} h y'(x_n) \\ &+ \underbrace{\left( \frac{1}{2} - \frac{1}{2}\alpha_1 + \beta_1 \right)}_{c_2} h^2 y''(x_n) \\ &+ \underbrace{\left( \frac{1}{6} + \frac{1}{6}\alpha_1 - \frac{1}{2}\beta_1 \right)}_{c_3} h^3 y'''(x_n) \\ &+ \underbrace{\left( \frac{1}{24} - \frac{1}{24}\alpha_1 + \frac{1}{6}\beta_1 \right)}_{c_4} h^4 y^{(4)}(x_n) + O(h^5) \end{align*} \]

步骤4：令低阶系数为0，建立方程组求解参数

要使方法的精度阶尽可能高，需让低阶项的系数全部为0，即\(c_0=c_1=c_2=c_3=0\)，得到线性方程组：

\[\begin{cases} \alpha_0 + \alpha_1 = 1 \\ -\alpha_1 + \beta_0 + \beta_1 = 1 \\ \alpha_1 - 2\beta_1 = 1 \\ -\alpha_1 + 3\beta_1 = 1 \end{cases} \]

解方程组：

  1. 联立后两个方程：\(\alpha_1 - 2\beta_1 = 1\)，\(-\alpha_1 + 3\beta_1 = 1\)，相加得\(\beta_1=2\)，代入得\(\alpha_1=5\)；
  2. 代入第一个方程，得\(\alpha_0=1 - \alpha_1 = -4\)；
  3. 代入第二个方程，得\(\beta_0=1 + \alpha_1 - \beta_1 = 1+5-2=4\)。

最终参数为：\(\boldsymbol{\alpha_0=-4,\ \alpha_1=5,\ \beta_0=4,\ \beta_1=2}\)。

步骤5：确定精度阶与局部截断误差

将参数代入\(c_4\)，计算得：

\[c_4 = \frac{1}{24} - \frac{5}{24} + \frac{2}{6} = \frac{1 - 5 + 8}{24} = \frac{4}{24} = \frac{1}{6} \neq 0 \]

因此，局部截断误差为：

\[T_{n+1} = \frac{1}{6} h^4 y^{(4)}(x_n) + O(h^5) \]

根据精度阶的定义，\(T_{n+1}=O(h^{p+1})\)，此处\(p+1=4\)，即\(p=3\)，因此该方法为三阶精度。

最终二步法公式

\[y_{n+1} = -4y_n + 5y_{n-1} + 2h\left( 2f_n + f_{n-1} \right), \quad n=0,1,2,\dots \]

三、例9.8 四阶线性多步法的参数证明

1. 题目重述

证明：存在\(\alpha\)的一个值，使线性多步法

\[y_{n+1} + \alpha(y_n - y_{n-1}) - y_{n-2} = \frac{1}{2}(3+\alpha)h\left( f_n + f_{n-1} \right), \quad n=2,3,\dots \]

为四阶方法。

2. 核心思路

根据精度阶的定义：若方法的局部截断误差\(T_{n+1}=O(h^5)\)，则方法为四阶精度。我们仍通过泰勒展开法，令\(h^0\)到\(h^4\)的系数全为0，求解对应的\(\alpha\)值。

3. 完整推导过程

步骤1：写出局部截断误差的定义

整理公式，写出局部截断误差（假设前序步均为精确值，\(f_n=y'(x_n)\)，\(f_{n-1}=y'(x_{n-1})\)）：

\[T_{n+1} = y(x_{n+1}) + \alpha\left[ y(x_n) - y(x_{n-1}) \right] - y(x_{n-2}) - \frac{1}{2}(3+\alpha)h\left[ y'(x_n) + y'(x_{n-1}) \right] \]

其中\(x_{n+1}=x_n+h\)，\(x_{n-1}=x_n-h\)，\(x_{n-2}=x_n-2h\)。

步骤2：泰勒展开所有项

将所有函数在\(x_n\)处做泰勒展开，展开到\(h^5\)项：

  1. \(y(x_n+h) = y + h y' + \frac{h^2}{2}y'' + \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} + \frac{h^5}{120}y^{(5)} + O(h^6)\)
  2. \(y(x_n-h) = y - h y' + \frac{h^2}{2}y'' - \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} - \frac{h^5}{120}y^{(5)} + O(h^6)\)
  3. \(y(x_n-2h) = y - 2h y' + \frac{(2h)^2}{2}y'' - \frac{(2h)^3}{6}y''' + \frac{(2h)^4}{24}y^{(4)} - \frac{(2h)^5}{120}y^{(5)} + O(h^6)\)
    \(= y - 2h y' + 2h^2 y'' - \frac{4h^3}{3}y''' + \frac{2h^4}{3}y^{(4)} - \frac{4h^5}{15}y^{(5)} + O(h^6)\)
  4. \(y'(x_n-h) = y' - h y'' + \frac{h^2}{2}y''' - \frac{h^3}{6}y^{(4)} + \frac{h^4}{24}y^{(5)} + O(h^5)\)

步骤3：代入并合并同类项

将展开式代入\(T_{n+1}\)，按\(h\)的幂次整理，合并同类项：

\[\begin{align*} T_{n+1} &= \left( y + h y' + \frac{h^2}{2}y'' + \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} \right) \\ &+ \alpha\left[ y - \left( y - h y' + \frac{h^2}{2}y'' - \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} \right) \right] \\ &- \left( y - 2h y' + 2h^2 y'' - \frac{4h^3}{3}y''' + \frac{2h^4}{3}y^{(4)} \right) \\ &- \frac{h}{2}(3+\alpha)\left[ y' + \left( y' - h y'' + \frac{h^2}{2}y''' - \frac{h^3}{6}y^{(4)} \right) \right] + O(h^5) \end{align*} \]

逐项化简，合并\(h\)的各阶系数：

  1. \(h^0\)项（\(y(x_n)\)的系数）：\(1 + \alpha - \alpha - 1 = 0\)，自动为0；
  2. \(h^1\)项（\(y'(x_n)\)的系数）：\(1 + \alpha + 2 - (3+\alpha) = 0\)，自动为0；
  3. \(h^2\)项（\(y''(x_n)\)的系数）：\(\frac{1}{2} - \frac{\alpha}{2} - 2 + \frac{3+\alpha}{2} = 0\)，自动为0；
  4. \(h^3\)项（\(y'''(x_n)\)的系数）：\(\frac{1}{6} + \frac{\alpha}{6} + \frac{4}{3} - \frac{3+\alpha}{4} = \frac{3}{4} - \frac{\alpha}{12}\)；
  5. \(h^4\)项（\(y^{(4)}(x_n)\)的系数）：\(\frac{1}{24} - \frac{\alpha}{24} - \frac{2}{3} + \frac{3+\alpha}{12} = \frac{1}{24}(-9 + \alpha)\)。

因此，局部截断误差可简化为：

\[T_{n+1} = \left( \frac{3}{4} - \frac{\alpha}{12} \right) h^3 y'''(x_n) + \frac{1}{24}(-9 + \alpha) h^4 y^{(4)}(x_n) + O(h^5) \]

步骤4：求解α使方法为四阶

要使方法为四阶精度，需满足\(T_{n+1}=O(h^5)\)，即\(h^3\)和\(h^4\)的系数全为0：

  1. 令\(h^3\)的系数为0：\(\frac{3}{4} - \frac{\alpha}{12} = 0 \implies \alpha=9\)；
  2. 代入\(\alpha=9\)，验证\(h^4\)的系数：\(\frac{1}{24}(-9+9)=0\)，恰好为0。

当\(\alpha=9\)时，\(T_{n+1}=O(h^5)\)，因此该方法为四阶精度，得证。

四、泰勒展开法构造多步法的通用步骤总结

  1. 定义局部截断误差：假设前序步数值解均为精确值，将\(f_{n+i}\)替换为\(y'(x_{n+i})\)，写出\(T_{n+k}\)的表达式；
  2. 泰勒展开：将所有\(y(x_n+ih)\)、\(y'(x_n+ih)\)在\(x_n\)处做泰勒展开，展开阶数需高于预期的精度阶；
  3. 合并同类项：按\(h\)的幂次整理\(T_{n+k}\)，得到各阶项的系数；
  4. 建立方程组求解参数：令低阶项的系数全为0，建立线性方程组，求解待定参数；
  5. 确定精度阶与误差：代入参数计算首个非零系数的阶数，确定精度阶，写出局部截断误差的主项。

该方法的核心优势是通用性极强，无论多步法的形式是否符合标准线性k步法，都可以通过该方法精准推导，是考试与工程应用的核心方法。

posted on 2026-03-07 23:33  Indian_Mysore  阅读(0)  评论(0)    收藏  举报

导航