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

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

 

ch01常微分方程初,边值问题数值解法

常微分方程初值问题数值解法 核心知识点详解与推导证明

我将以多年教学与研究的经验，从背景与问题定义、核心理论定理（含完整证明）、数值方法基本思想、方法分类与特性四个维度展开讲解，最终用表格完成系统性归纳总结。

一、引言与问题背景

常微分方程（ODE）与微积分同时诞生，Newton、Leibniz、Euler等先驱早在17-18世纪就开展了相关研究，但工程与科学实践中，绝大多数常微分方程无法通过初等积分法得到解析解，仅有极少数方程能写出显式通解。

这一局限性催生了两个研究方向：一是从理论上研究解的定性性质，二是从工程角度研究解的近似计算方法。随着计算机技术的发展，数值解法成为求解常微分方程最核心、最实用的手段，它的核心目标是：不求解析表达式，直接计算解函数在离散节点上的近似值。

我们的研究核心是一阶常微分方程初值问题（IVP），其标准形式为：

\[\begin{cases} \frac{dy}{dx} = f(x,y), \quad x \in [a,b] \tag{1.1.1} \\ y(a) = y_a \tag{1.1.2} \end{cases} \]

其中：

  • \(f(x,y)\)是定义在区域\(\Omega = \{(x,y)\mid x\in [a,b],\ y\in \mathbb{R}\}\)上的二元实值函数，描述解的变化率；
  • \(y_a\)是给定的初始条件，即解在区间左端点\(x=a\)处的函数值；
  • 我们需要求解的是：在\([a,b]\)上一阶连续可微的函数\(y(x)\)，使其同时满足微分方程与初始条件。

二、初值问题的核心理论定理与完整证明

数值解法的前提是：问题本身是“数学上合理的”——即解存在、唯一，且对数据扰动不敏感。这部分我们将完整证明两个核心定理：解的存在唯一性定理与适定性定理。

2.1 解的存在唯一性定理（定理1.1.1）

定理内容

若方程(1.1.1)的右端函数\(f(x,y)\)满足以下3个条件：

  1. \(f(x,y)\)是实值函数；
  2. \(f(x,y)\)在矩形区域\(\Omega = [a,b] \times \mathbb{R}\)上连续；
  3. \(f(x,y)\)关于\(y\)满足Lipschitz条件：存在正常数\(L\)（Lipschitz常数），使得对任意\(x\in [a,b]\)、任意\(y,z\in \mathbb{R}\)，均有

    \[|f(x,y) - f(x,z)| \leq L |y - z| \tag{1.1.3} \]

则初值问题(1.1.1)-(1.1.2)存在唯一的解\(y(x) \in C^1[a,b]\)（即\(y(x)\)在\([a,b]\)上一阶连续可微）。

条件解读

  1. 连续性条件：是解存在的基础（Peano存在定理指出，仅需\(f\)连续，解就至少存在一个），但无法保证唯一性；
  2. Lipschitz条件：是解唯一性的核心，它是比“关于\(y\)可微”更弱、比“关于\(y\)连续”更强的条件：
    • 若\(\frac{\partial f}{\partial y}\)在\(\Omega\)上存在且有界（\(|\frac{\partial f}{\partial y}| \leq L\)），由微分中值定理可直接推出Lipschitz条件；
    • 满足Lipschitz条件的函数不一定可微（例如\(f(x,y)=|y|\)，Lipschitz常数\(L=1\)，但在\(y=0\)处不可微）。

完整证明（Picard迭代法，构造性证明）

我们将证明分为5个核心步骤，该方法同时也是数值迭代思想的源头。

步骤1：将初值问题转化为等价积分方程

若\(y(x)\)是初值问题的解，则\(y'(x)=f(x,y(x))\)，对等式两边从\(a\)到\(x\)积分：

\[\int_{a}^x y'(t) dt = \int_{a}^x f(t,y(t)) dt \]

由牛顿-莱布尼茨公式，左边为\(y(x)-y(a)\)，代入初始条件\(y(a)=y_a\)，得：

\[y(x) = y_a + \int_{a}^x f(t,y(t)) dt \tag{1} \]

反之，若\(y(x)\)满足积分方程(1)，两边对\(x\)求导（变上限积分求导法则）得\(y'(x)=f(x,y(x))\)，代入\(x=a\)得\(y(a)=y_a\)，即满足初值问题。

因此，初值问题与积分方程(1)完全等价，我们只需证明积分方程(1)存在唯一连续解。

步骤2：构造Picard迭代序列

我们构造逐步逼近的函数序列\(\{y_n(x)\}\)，迭代格式为：

\[\begin{cases} y_0(x) = y_a, \quad \forall x\in [a,b] \quad \text{（零次迭代，初始常数函数）} \\ y_{n+1}(x) = y_a + \int_{a}^x f(t,y_n(t)) dt, \quad n=0,1,2,\dots \end{cases} \]

我们将证明该序列一致收敛，且极限就是积分方程的解。

步骤3：证明迭代序列在\([a,b]\)上一致收敛

将序列改写为级数形式：

\[y_n(x) = y_0(x) + \sum_{k=0}^{n-1} \left[ y_{k+1}(x) - y_k(x) \right] \]

序列的一致收敛性等价于上述函数项级数的一致收敛性，我们用Weierstrass优级数判别法证明。

首先用数学归纳法证明通项的上界：

  • 当\(k=0\)时：

    \[|y_1(x)-y_0(x)| = \left| \int_{a}^x f(t,y_a) dt \right| \]

    闭区域上的连续函数必有界，设\(M = \max_{(x,y)\in \Omega} |f(x,y)|\)，因此：

    \[|y_1(x)-y_0(x)| \leq \int_{a}^x |f(t,y_a)| dt \leq M (x-a) \leq M(b-a) \]

  • 归纳假设：当\(k=m\)时，有

    \[|y_m(x) - y_{m-1}(x)| \leq \frac{M L^{m-1} (x-a)^m}{m!}, \quad \forall x\in [a,b] \]

  • 归纳递推：当\(k=m+1\)时：

    \[|y_{m+1}(x)-y_m(x)| = \left| \int_{a}^x \left[ f(t,y_m(t)) - f(t,y_{m-1}(t)) \right] dt \right| \]

    由Lipschitz条件，\(|f(t,y_m)-f(t,y_{m-1})| \leq L |y_m(t)-y_{m-1}(t)|\)，代入归纳假设：

    \[|y_{m+1}(x)-y_m(x)| \leq L \int_{a}^x \frac{M L^{m-1} (t-a)^m}{m!} dt \]

    计算积分得：

    \[|y_{m+1}(x)-y_m(x)| \leq \frac{M L^m (x-a)^{m+1}}{(m+1)!} \]

由数学归纳法，对所有\(n\geq0\)，均有：

\[|y_{n+1}(x)-y_n(x)| \leq \frac{M L^n (b-a)^{n+1}}{(n+1)!}, \quad \forall x\in [a,b] \]

考察优级数\(\sum_{n=0}^\infty \frac{M L^n (b-a)^{n+1}}{(n+1)!}\)，其和为\(\frac{M}{L}\left( e^{L(b-a)} - 1 \right)\)，是收敛的正项级数。由Weierstrass判别法，函数项级数在\([a,b]\)上一致收敛，因此迭代序列\(\{y_n(x)\}\)在\([a,b]\)上一致收敛，记极限函数为：

\[\lim_{n\to\infty} y_n(x) = y(x), \quad x\in [a,b] \]

步骤4：证明极限函数是积分方程的解
  1. 连续性：连续函数序列一致收敛的极限仍为连续函数，因此\(y(x)\)在\([a,b]\)上连续；
  2. 满足积分方程：对迭代式两边取\(n\to\infty\)的极限：

    \[\lim_{n\to\infty} y_{n+1}(x) = y_a + \lim_{n\to\infty} \int_{a}^x f(t,y_n(t)) dt \]

    由\(f\)的连续性与\(y_n(t)\)的一致收敛性，\(f(t,y_n(t))\)一致收敛到\(f(t,y(t))\)，因此极限与积分可交换顺序：

    \[\lim_{n\to\infty} \int_{a}^x f(t,y_n(t)) dt = \int_{a}^x \lim_{n\to\infty} f(t,y_n(t)) dt = \int_{a}^x f(t,y(t)) dt \]

    因此得到：

    \[y(x) = y_a + \int_{a}^x f(t,y(t)) dt \]

    即\(y(x)\)是积分方程的解，也是初值问题的解，且\(y(x)\in C^1[a,b]\)。
步骤5：证明解的唯一性

假设积分方程有两个连续解\(y(x)\)和\(z(x)\)，则：

\[y(x)-z(x) = \int_{a}^x \left[ f(t,y(t)) - f(t,z(t)) \right] dt \]

两边取绝对值，结合Lipschitz条件得：

\[|y(x)-z(x)| \leq L \int_{a}^x |y(t)-z(t)| dt \tag{2} \]

设\(K = \max_{x\in [a,b]} |y(x)-z(x)|\)，代入(2)得\(|y(x)-z(x)| \leq LK(x-a)\)，将该结果再次代入(2)，反复迭代\(n\)次后得：

\[|y(x)-z(x)| \leq K \cdot \frac{[L(b-a)]^n}{n!} \]

由于\(\lim_{n\to\infty} \frac{[L(b-a)]^n}{n!} = 0\)，令\(n\to\infty\)得\(|y(x)-z(x)| \leq 0\)，因此\(y(x)=z(x)\)，解唯一。

至此，存在唯一性定理全部证明完毕。

2.2 初值问题的适定性

适定性的定义（定义1.1.1）

称初值问题(1.1.1)-(1.1.2)对初值\(y_a\)是适定的，如果存在常数\(K>0\)和\(\eta>0\)，使得对任意\(0<\varepsilon\leq\eta\)，当

\[|y_a - \tilde{y}_a| < \varepsilon, \quad |f(x,y) - \tilde{f}(x,y)| < \varepsilon, \quad \forall (x,y)\in \Omega \]

时，扰动后的初值问题

\[\begin{cases} \frac{dz}{dx} = \tilde{f}(x,z), \quad x\in [a,b] \\ z(a) = \tilde{y}_a \end{cases} \]

存在解\(z(x)\)，且满足\(|y(x) - z(x)| \leq K\varepsilon\)，其中\(y(x)\)是原问题的解。

定义解读

适定性是Hadamard提出的核心概念，一个适定的问题需满足三个要求：

  1. 解存在；2. 解唯一；3. 解连续依赖于初始数据和方程右端项。
    数值计算中必然存在舍入误差、截断误差，这些都可看作对原问题的微小扰动。若问题不适定，微小扰动会导致解的巨大偏差，数值计算将失去意义。

适定性定理（定理1.1.2）

若方程右端函数\(f(x,y)\)在区域\(\Omega\)上满足Lipschitz条件(1.1.3)，则初值问题(1.1.1)-(1.1.2)对任何初值\(y_a\)都是适定的。

完整证明

原问题的解\(y(x)\)与扰动问题的解\(z(x)\)分别满足积分方程：

\[y(x) = y_a + \int_{a}^x f(t,y(t)) dt, \quad z(x) = \tilde{y}_a + \int_{a}^x \tilde{f}(t,z(t)) dt \]

两式相减取绝对值，由三角不等式得：

\[|y(x)-z(x)| \leq |y_a - \tilde{y}_a| + \left| \int_{a}^x \left[ f(t,y(t)) - \tilde{f}(t,z(t)) \right] dt \right| \]

将被积函数拆分为两项：

\[f(t,y(t)) - \tilde{f}(t,z(t)) = \left[ f(t,y(t)) - f(t,z(t)) \right] + \left[ f(t,z(t)) - \tilde{f}(t,z(t)) \right] \]

结合Lipschitz条件\(|f(t,y)-f(t,z)| \leq L|y-z|\)，以及扰动条件\(|y_a-\tilde{y}_a|<\varepsilon\)、\(|f-\tilde{f}|<\varepsilon\)，代入得：

\[|y(x)-z(x)| < \varepsilon + \int_{a}^x L|y(t)-z(t)| dt + \int_{a}^x \varepsilon dt \]

计算积分项\(\int_{a}^x \varepsilon dt = \varepsilon(x-a) \leq \varepsilon(b-a)\)，整理得：

\[|y(x)-z(x)| < \varepsilon\left[ 1 + (b-a) \right] + L \int_{a}^x |y(t)-z(t)| dt \tag{3} \]

此时我们引入Gronwall积分不等式（先证明该不等式，再代入求解）：

Gronwall不等式：设\(u(x),k(x)\)是\([a,b]\)上的非负连续函数，\(C\geq0\)为常数，且满足

\[u(x) \leq C + \int_{a}^x k(t)u(t) dt, \quad \forall x\in [a,b] \]

则有

\[u(x) \leq C \exp\left( \int_{a}^x k(t) dt \right), \quad \forall x\in [a,b] \]

Gronwall不等式证明：
令\(U(x) = C + \int_{a}^x k(t)u(t) dt\)，则\(u(x)\leq U(x)\)，且\(U'(x)=k(x)u(x)\leq k(x)U(x)\)，\(U(a)=C\)。
当\(C>0\)时，对\(\ln U(x)\)求导得：

\[\frac{d}{dx}\ln U(x) = \frac{U'(x)}{U(x)} \leq k(x) \]

两边从\(a\)到\(x\)积分得：

\[\ln U(x) - \ln C \leq \int_{a}^x k(t) dt \implies U(x) \leq C \exp\left( \int_{a}^x k(t) dt \right) \]

结合\(u(x)\leq U(x)\)，不等式得证。\(C=0\)时可通过极限方法同理得证。

回到不等式(3)，令\(u(x)=|y(x)-z(x)|\)，\(C=\varepsilon\left[1+(b-a)\right]\)，\(k(t)=L\)，代入Gronwall不等式得：

\[|y(x)-z(x)| < \varepsilon\left[1+(b-a)\right] \cdot e^{L(x-a)} \]

由于\(x\in [a,b]\)，\(e^{L(x-a)} \leq e^{L(b-a)}\)，因此：

\[|y(x)-z(x)| < \varepsilon \cdot \underbrace{\left[1+(b-a)\right]e^{L(b-a)}}_{K} \]

令\(K=\left[1+(b-a)\right]e^{L(b-a)}\)，显然\(K>0\)且与\(\varepsilon\)无关，满足适定性定义的要求，因此初值问题是适定的。

三、数值解法的基本思想与分类

3.1 数值解法的核心：离散化

初值问题的解\(y(x)\)是连续区间\([a,b]\)上的函数，数值解法的核心是离散化：

  1. 对区间\([a,b]\)做网格剖分，取离散节点：\(a = x_0 < x_1 < x_2 < \dots < x_N = b\)；
  2. 通常取等间距节点，步长\(h = \frac{b-a}{N}\)，此时\(x_m = a + mh\)，\(m=0,1,\dots,N\)；
  3. 目标是计算解在节点处的近似值\(y_m \approx y(x_m)\)，\(m=1,2,\dots,N\)（\(y_0=y_a\)已知）；
  4. 将连续的微分方程转化为离散节点上的代数方程（差分方程），通过递推求解近似值。

常用离散化方法

  1. 差商近似微商法：用有限差商代替导数（如向前差商\(\frac{y(x_{m+1})-y(x_m)}{h} \approx y'(x_m)\)），直接将微分方程转化为差分方程；
  2. Taylor级数展开法：将\(y(x_{m+1})\)在\(x_m\)处做Taylor展开，截断到指定阶数，构造高阶精度的递推公式；
  3. 数值积分法：将微分方程在\([x_m,x_{m+1}]\)上积分，用数值积分公式近似积分项，得到离散递推关系。

3.2 数值方法的核心分类

根据计算\(y_{m+1}\)时用到的节点数量，数值方法分为两大类：单步法与多步法。

1. 单步法

定义：计算\(x_{m+1}\)处的近似值\(y_{m+1}\)时，仅用到前一个节点\(x_m\)处的信息（\(y_m\)、\(f(x_m,y_m)\)等），无需更早节点的信息。
核心特点：

  • 可自启动：仅需初始值\(y_0\)，即可依次计算\(y_1,y_2,\dots,y_N\)；
  • 步长调整灵活：计算过程中可随时改变步长，无需额外处理；
  • 高阶单步法每一步需要多次计算\(f(x,y)\)的值。
    典型方法：Euler方法（显式/隐式）、梯形法、Runge-Kutta（龙格-库塔）系列方法。

2. 多步法

定义：计算\(y_{m+1}\)时，需要用到\(x_{m+1}\)左侧的多个节点（如\(x_m,x_{m-1},\dots,x_{m-k+1}\)，共\(k\)个节点，称为\(k\)步法）处的\(y\)值与\(f\)值。
核心特点：

  • 不可自启动：\(k\)步法需要先通过单步法计算前\(k\)个节点的近似值\(y_0,y_1,\dots,y_{k-1}\)，才能启动递推；
  • 步长调整较复杂：改变步长时需要重新计算初始节点值或做插值处理；
  • 同等精度下，计算量远小于单步法：线性多步法每一步仅需计算1次\(f(x,y)\)的值。
    典型方法：Adams显式/隐式方法、Milne方法、Hamming方法。

补充：显式方法与隐式方法

  • 显式方法：\(y_{m+1}\)可通过已知量直接计算，无需求解方程，计算简单、速度快，但稳定性相对较弱；
  • 隐式方法：\(y_{m+1}\)出现在等式两端（如\(f\)的自变量中），需要迭代求解方程，稳定性好、数值阻尼强，适合求解刚性方程。

四、核心知识点系统性归纳总结

表1 初值问题核心定理总结

定理名称 核心条件 核心结论 核心意义
存在唯一性定理 1. \(f(x,y)\)是\(\Omega=[a,b]\times\mathbb{R}\)上的实值连续函数；
2. \(f(x,y)\)关于\(y\)满足Lipschitz条件		 初值问题存在唯一的一阶连续可微解\(y(x)\in C^1[a,b]\) 1. 为数值求解提供了理论前提，保证解的存在性与唯一性；
2. 构造性的Picard迭代为数值迭代方法提供了核心思想
适定性定理 \(f(x,y)\)在\(\Omega\)上满足Lipschitz条件 初值问题对任意初值都是适定的，解连续依赖于初始值与方程右端项，微小扰动仅引起解的微小变化 1. 保证数值计算的可靠性，避免误差被无限放大；
2. 是数值方法收敛性、稳定性分析的核心基础

表2 离散化方法总结

离散化方法 核心思想 典型应用 核心特点
差商近似微商法 用有限差商代替微分方程中的导数，直接转化为差分方程 Euler方法、向后Euler方法 直观易懂，推导简单，是最基础的离散化方法
Taylor级数展开法 对解做Taylor展开并截断到指定阶数，构造递推公式 Runge-Kutta方法、Taylor级数法 可构造任意高阶精度的方法，是高阶单步法的核心构造手段
数值积分法 对微分方程做区间积分，用数值积分公式近似积分项 Adams多步法、梯形法 可直接通过数值积分的精度分析得到方法的阶数，是线性多步法的核心构造手段

表3 数值方法分类与特性总结

方法类别 核心定义 核心特点 典型方法 适用场景
单步法 计算\(y_{m+1}\)仅用到前一个节点\(x_m\)的信息 1. 可自启动，仅需初始值\(y_0\)；
2. 步长调整灵活；
3. 高阶方法每步需多次计算\(f\)值		 显式/隐式Euler法、梯形法、Runge-Kutta系列方法 步长需频繁调整的问题、非刚性方程、需要自启动的场景
多步法 计算\(y_{m+1}\)需要用到多个历史节点的信息 1. 不可自启动，需单步法提供初始值；
2. 步长调整较复杂；
3. 同等精度下计算量更小		 Adams显式/隐式方法、Milne方法、Hamming方法 大规模长区间求解、对计算效率要求高的场景
显式方法 \(y_{m+1}\)可由已知量直接计算，无需求解方程 计算简单、速度快；稳定性相对较弱 显式Euler法、Adams-Bashforth方法、显式RK方法 非刚性方程、步长较小的常规求解场景
隐式方法 \(y_{m+1}\)出现在方程两端，需迭代求解 稳定性好、数值阻尼强；计算相对复杂 隐式Euler法、梯形法、Adams-Moulton方法、隐式RK方法 刚性方程、对稳定性要求高的场景、大步长计算

Euler方法 完整讲解与推导证明

承接上一章的一阶常微分方程初值问题

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

我们从最基础、最经典的数值解法——Euler方法（欧拉方法） 展开讲解。Euler方法是所有常微分方程数值解法的基石，其核心思想是“以直代曲”，通过离散化将连续微分方程转化为可递推的代数差分方程，兼具直观的几何意义与严谨的理论基础。

一、Euler方法的基本格式与几何意义

1.1 基本格式的推导

我们从两个等价的角度推导Euler格式，帮助理解其本质：

角度1：向前差商近似微商

导数的定义是：

\[y'(x_m) = \lim_{h \to 0} \frac{y(x_{m+1}) - y(x_m)}{h} \]

当步长\(h\)足够小时，我们用向前差商近似代替导数，即：

\[\frac{y(x_{m+1}) - y(x_m)}{h} \approx y'(x_m) = f(x_m, y(x_m)) \]

将上式整理，用近似值\(y_m\)代替精确值\(y(x_m)\)，\(y_{m+1}\)代替\(y(x_{m+1})\)，得到递推公式：

\[\boldsymbol{y_{m+1} = y_m + h f(x_m, y_m), \quad m=0,1,\dots,N-1} \tag{1.2.1} \]

这就是显式Euler方法的核心递推格式。其中\(h = \frac{b-a}{N}\)为等距步长，\(x_m = a + mh\)为离散节点。

角度2：数值积分的矩形公式近似

对微分方程\(y'(s) = f(s,y(s))\)在区间\([x_m, x_{m+1}]\)上积分，得到精确的积分等式：

\[y(x_{m+1}) = y(x_m) + \int_{x_m}^{x_{m+1}} f(s,y(s)) ds \tag{1.2.2} \]

这是所有单步法的核心等式，不同的数值积分方法对应不同的单步格式。

对于Euler方法，我们采用左矩形数值积分公式近似右端的积分项：

\[\int_{x_m}^{x_{m+1}} f(s,y(s)) ds \approx h \cdot f(x_m, y(x_m)) \]

代入积分等式，同样得到Euler格式(1.2.1)。

1.2 几何意义（折线法）

Euler方法有非常直观的几何解释，因此也被称为折线法：

  1. 初值问题的精确解\(y(x)\)是\((x,y)\)平面上过初始点\((x_0,y_0)\)的一条光滑积分曲线，该曲线上任意一点\((x,y)\)的切线斜率为\(f(x,y)\)。
  2. 在初始点\((x_0,y_0)\)处，我们作积分曲线的切线，其斜率为\(f(x_0,y_0)\)，切线方程为\(y = y_0 + f(x_0,y_0)(x-x_0)\)。
  3. 取\(x=x_1=x_0+h\)，切线上的点的纵坐标为\(y_1 = y_0 + h f(x_0,y_0)\)，我们用\(y_1\)作为精确值\(y(x_1)\)的近似值。
  4. 接着，在点\((x_1,y_1)\)处，再次作斜率为\(f(x_1,y_1)\)的切线，取\(x=x_2=x_1+h\)，得到\(y_2 = y_1 + h f(x_1,y_1)\)，作为\(y(x_2)\)的近似值。
  5. 重复上述过程，最终得到一条折线\((x_0,y_0) \to (x_1,y_1) \to \dots \to (x_N,y_N)\)，我们用这条折线来近似代替原光滑的积分曲线。

这就是Euler方法的几何本质：用分段的切线折线，逐段近似光滑的积分曲线。

二、Euler方法的误差分析

数值方法的核心性能指标之一是精度，我们通过截断误差来刻画Euler方法的精度，分为局部截断误差和整体截断误差。

2.1 核心概念定义

首先明确两个关键误差的定义：

  1. 局部截断误差\(R_m\)：假设前一步的计算值是精确的（即\(y_m = y(x_m)\)），用数值方法计算一步得到的\(y_{m+1}\)与精确值\(y(x_{m+1})\)的误差，即：

    \[R_m = y(x_{m+1}) - \left[ y(x_m) + h f(x_m, y(x_m)) \right] \]

    它刻画的是单步计算本身带来的误差，不包含之前步骤的误差累积。

  2. 整体截断误差\(e_m\)：考虑所有前序步骤的误差累积，计算值\(y_m\)与精确值\(y(x_m)\)的总误差，即：

    \[e_m = y(x_m) - y_m \]

    它刻画的是递推到第\(m\)步时的总误差，是我们最终关心的实际计算误差。

2.2 局部截断误差的估计（引理1.2.1）

引理内容

设\(f(x,y)\)关于\(x\)和\(y\)均满足Lipschitz条件，对应的Lipschitz常数分别为\(K\)和\(L\)，记\(M = \max_{x\in [a,b]} |y'(x)| = \max_{x\in [a,b]} |f(x,y(x))|\)，则Euler方法的局部截断误差\(R_m\)满足：

\[|R_m| \leq \frac{h^2}{2}(K + LM), \quad x_m < \xi < x_{m+1} \tag{1.2.5} \]

即Euler方法的局部截断误差是\(O(h^2)\)阶的。

完整证明

从局部截断误差的定义出发：

\[R_m = y(x_{m+1}) - y(x_m) - h f(x_m, y(x_m)) \]

首先，\(y(x_{m+1}) - y(x_m) = \int_{x_m}^{x_{m+1}} y'(s) ds = \int_{x_m}^{x_{m+1}} f(s,y(s)) ds\)，因此：

\[R_m = \int_{x_m}^{x_{m+1}} f(s,y(s)) ds - h f(x_m, y(x_m)) \tag{1.2.4} \]

为了估计这个积分的上界，我们对被积函数做拆分，利用三角不等式：

\[\begin{align*} |R_m| &= \left| \int_{x_m}^{x_{m+1}} \left[ f(s,y(s)) - f(x_m, y(x_m)) \right] ds \right| \\ &= \left| \int_{x_m}^{x_{m+1}} \left[ f(s,y(s)) - f(x_m, y(s)) + f(x_m, y(s)) - f(x_m, y(x_m)) \right] ds \right| \\ &\leq \int_{x_m}^{x_{m+1}} \left| f(s,y(s)) - f(x_m, y(s)) \right| ds + \int_{x_m}^{x_{m+1}} \left| f(x_m, y(s)) - f(x_m, y(x_m)) \right| ds \end{align*} \]

接下来分别估计两个积分项：

  1. 第一项估计：利用\(f\)关于\(x\)的Lipschitz条件，\(|f(s,y) - f(x_m,y)| \leq K |s - x_m|\)，因此：

    \[\int_{x_m}^{x_{m+1}} \left| f(s,y(s)) - f(x_m, y(s)) \right| ds \leq K \int_{x_m}^{x_{m+1}} |s - x_m| ds \]

    计算积分：\(\int_{x_m}^{x_{m+1}} (s - x_m) ds = \frac{1}{2}h^2\)，因此第一项上界为\(\frac{1}{2} K h^2\)。

  2. 第二项估计：利用\(f\)关于\(y\)的Lipschitz条件，\(|f(x_m,y(s)) - f(x_m,y(x_m))| \leq L |y(s) - y(x_m)|\)。
    而\(y(s) - y(x_m) = \int_{x_m}^s y'(t) dt\)，因此\(|y(s)-y(x_m)| \leq \int_{x_m}^s |y'(t)| dt \leq M (s - x_m)\)，代入得：

    \[\int_{x_m}^{x_{m+1}} \left| f(x_m, y(s)) - f(x_m, y(x_m)) \right| ds \leq L M \int_{x_m}^{x_{m+1}} (s - x_m) ds = \frac{1}{2} L M h^2 \]

将两项的上界合并，得到：

\[|R_m| \leq \frac{1}{2} K h^2 + \frac{1}{2} L M h^2 = \frac{h^2}{2}(K + LM) \]

引理得证。

补充直观推导：若\(y(x)\)二阶连续可微，将\(y(x_{m+1})\)在\(x_m\)处做Taylor展开：

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

代入\(y'(x_m)=f(x_m,y(x_m))\)，与Euler格式对比，直接得到\(R_m = \frac{h^2}{2} y''(\xi) = O(h^2)\)，与引理结论完全一致。

2.3 整体截断误差的估计（定理1.2.1）

定理内容

设\(f(x,y)\)关于\(x\)和\(y\)均满足Lipschitz条件，\(K\)和\(L\)为对应的Lipschitz常数，且当\(h \to 0\)时，\(y_0 \to y(x_0)\)，则Euler方法的解\(\{y_m\}\)一致收敛于初值问题的精确解，且整体截断误差\(e_m\)满足估计式：

\[|e_m| \leq e^{L(b-a)} |e_0| + \frac{h}{2} \left( \frac{K}{L} + M \right) \left( e^{L(b-a)} - 1 \right) \tag{1.2.6} \]

若初始值精确，即\(e_0 = y(x_0) - y_0 = 0\)，则：

\[|e_m| = O(h) \]

即Euler方法的整体截断误差是\(O(h)\)阶的，比局部截断误差低一阶。

完整证明

首先推导整体截断误差的递推关系。
精确解满足：

\[y(x_{m+1}) = y(x_m) + h f(x_m, y(x_m)) + R_m \]

Euler方法的计算值满足：

\[y_{m+1} = y_m + h f(x_m, y_m) \]

两式相减，得到整体截断误差的递推式：

\[e_{m+1} = e_m + h \left[ f(x_m, y(x_m)) - f(x_m, y_m) \right] + R_m \tag{1.2.7} \]

对等式两边取绝对值，利用三角不等式和\(f\)关于\(y\)的Lipschitz条件\(|f(x_m,y(x_m)) - f(x_m,y_m)| \leq L |e_m|\)，得：

\[|e_{m+1}| \leq |e_m| + h L |e_m| + |R_m| = (1 + hL) |e_m| + |R_m| \]

代入局部截断误差的上界\(|R_m| \leq R = \frac{h^2}{2}(K + LM)\)，得到：

\[|e_{m+1}| \leq (1 + hL) |e_m| + R \]

对该线性递推不等式展开递推：

  • \(k=1\)时：\(|e_1| \leq (1+hL)|e_0| + R\)
  • \(k=2\)时：\(|e_2| \leq (1+hL)|e_1| + R \leq (1+hL)^2 |e_0| + R \left[ 1 + (1+hL) \right]\)
  • 递推到第\(k\)步时：

    \[|e_k| \leq (1+hL)^k |e_0| + R \sum_{i=0}^{k-1} (1+hL)^i \]

右边求和项为等比数列，公比\(1+hL\)，求和得：

\[\sum_{i=0}^{k-1} (1+hL)^i = \frac{(1+hL)^k - 1}{hL} \]

代入递推式，结合不等式\(1 + t \leq e^t\)（\(t>0\)），令\(t=hL\)，得\((1+hL)^k \leq e^{khL} \leq e^{L(b-a)}\)（因\(kh \leq b-a\)）。

将\(R = \frac{h^2}{2}(K + LM)\)代入化简\(\frac{R}{hL}\)：

\[\frac{R}{hL} = \frac{h}{2} \left( \frac{K}{L} + M \right) \]

最终得到整体截断误差的上界：

\[|e_k| \leq e^{L(b-a)} |e_0| + \frac{h}{2} \left( \frac{K}{L} + M \right) \left( e^{L(b-a)} - 1 \right) \]

定理得证。

关键结论：初始值精确时，\(|e_m|\)与\(h\)成正比，整体截断误差为一阶。这是因为每一步的局部二阶误差会随步数累积，最终总误差阶数降低一阶；当\(h \to 0\)时，\(|e_m| \to 0\)，说明Euler方法是收敛的。

三、Euler方法的稳定性

数值方法的稳定性，刻画的是初始值扰动、舍入误差等，是否会随计算步数无限放大。若误差不会被放大，方法是稳定的，否则为不稳定。

3.1 稳定性的定义（定义1.2.1）

称Euler方法是稳定的（零稳定性，关于初值的稳定性），如果存在常数\(C>0\)和\(h_0>0\)，使得对任意\(0 < h < h_0\)、任意\(a \leq mh \leq b\)，Euler方法的两个解\(y_m\)和\(z_m\)满足：

\[|y_m - z_m| \leq C |y_0 - z_0| \tag{1.2.8} \]

其中\(y_m\)和\(z_m\)分别是以\(y_0\)和\(z_0\)为初始值的Euler方法的解。

核心含义：初始值的微小扰动，在整个计算过程中不会被无限放大，始终被初始扰动的常数倍控制。

3.2 稳定性定理（定理1.2.3）

定理内容

在\(f(x,y)\)关于\(y\)满足Lipschitz条件的假设下，Euler方法是稳定的。

完整证明

设\(y_m\)和\(z_m\)分别是初始值为\(y_0\)和\(z_0\)的Euler解，满足：

\[y_{m+1} = y_m + h f(x_m, y_m) \]

\[z_{m+1} = z_m + h f(x_m, z_m) \]

两式相减，记\(e_m = y_m - z_m\)，得到误差递推式：

\[e_{m+1} = e_m + h \left[ f(x_m, y_m) - f(x_m, z_m) \right] \]

两边取绝对值，结合\(f\)的Lipschitz条件\(|f(x_m,y_m)-f(x_m,z_m)| \leq L |e_m|\)，得：

\[|e_{m+1}| \leq (1 + hL) |e_m| \]

递推展开得：

\[|e_m| \leq (1 + hL)^m |e_0| \]

利用\(1+hL \leq e^{hL}\)，且\(mh \leq b-a\)，因此\((1 + hL)^m \leq e^{L(b-a)}\)。

令\(C = e^{L(b-a)}\)（与\(h\)、\(m\)无关的正常数），则：

\[|e_m| \leq C |e_0| \]

完全符合稳定性定义，定理得证。

四、改进的Euler方法

显式Euler方法结构简单但精度较低（整体一阶），我们通过更换更高精度的数值积分公式，得到改进格式。

4.1 梯形方法（隐式改进格式）

格式推导

回到核心积分等式，改用二阶梯形数值积分公式近似积分项：

\[\int_{x_m}^{x_{m+1}} f(s,y(s)) ds \approx \frac{h}{2} \left[ f(x_m, y(x_m)) + f(x_{m+1}, y(x_{m+1})) \right] \]

代入积分等式，得到梯形方法递推格式：

\[\boldsymbol{y_{m+1} = y_m + \frac{h}{2} \left[ f(x_m, y_m) + f(x_{m+1}, y_{m+1}) \right]} \tag{1.2.9} \]

格式特性

  1. 隐式格式：等式两端均含\(y_{m+1}\)，无法直接计算，需迭代求解；
  2. 精度：通过Taylor展开可证明，局部截断误差\(R_m^{(1)} = -\frac{h^3}{12} y'''(\xi) = O(h^3)\)，整体截断误差为\(O(h^2)\)，精度比显式Euler高一阶；
  3. 稳定性：绝对稳定区间为\((-\infty,0)\)，A-稳定，稳定性远优于显式Euler。

迭代求解

迭代格式为：

  1. 预测步：显式Euler提供初始值\(y_{m+1}^{(0)} = y_m + h f(x_m, y_m)\)
  2. 校正步：迭代更新\(y_{m+1}^{(k+1)} = y_m + \frac{h}{2} \left[ f(x_m, y_m) + f(x_{m+1}, y_{m+1}^{(k)}) \right]\)

当步长满足\(\frac{hL}{2} < 1\)时迭代收敛，通常迭代2~3次即可达到精度要求。

4.2 改进的Euler方法（预测-校正格式）

为避免隐式迭代，将预测步和校正步结合，仅做一次校正，得到工程中最常用的改进Euler方法：

\[\begin{cases} \text{预测步：} \quad \tilde{y}_{m+1} = y_m + h f(x_m, y_m) \\ \text{校正步：} \quad y_{m+1} = y_m + \frac{h}{2} \left[ f(x_m, y_m) + f(x_{m+1}, \tilde{y}_{m+1}) \right] \end{cases} \tag{1.2.12} \]

格式特性

  1. 显式格式：无需迭代，每步仅计算2次\(f(x,y)\)，计算效率高；
  2. 精度：整体截断误差\(O(h^2)\)，与梯形方法一致；
  3. 易用性：兼顾精度与计算效率，是入门级最实用的ODE数值解法。

五、核心知识点系统性归纳总结

表1 Euler系列方法核心格式与误差特性

方法名称 递推格式 局部截断误差 整体截断误差 格式类型 每步\(f\)计算次数
显式Euler方法 \(y_{m+1} = y_m + h f(x_m,y_m)\) \(O(h^2)\) \(O(h)\) 显式 1次
梯形方法 \(y_{m+1} = y_m + \frac{h}{2}\left[ f(x_m,y_m) + f(x_{m+1},y_{m+1}) \right]\) \(O(h^3)\) \(O(h^2)\) 隐式 迭代次数决定
改进Euler方法（预测-校正） 预测：\(\tilde{y}_{m+1}=y_m+hf(x_m,y_m)\)
校正：\(y_{m+1}=y_m+\frac{h}{2}\left[ f(x_m,y_m)+f(x_{m+1},\tilde{y}_{m+1}) \right]\)		 \(O(h^3)\) \(O(h^2)\) 显式 2次

表2 Euler系列方法稳定性与适用场景

方法名称 零稳定性 绝对稳定区间 核心优点 核心缺点 适用场景
显式Euler方法 稳定 \((-2,0)\) 结构简单、计算量极小、可自启动、步长调整灵活 精度低、稳定性较差 教学演示、步长极小的简单问题、非刚性方程快速求解
梯形方法 稳定 \((-\infty,0)\)（A-稳定） 精度高、稳定性极好、数值阻尼强 隐式格式需迭代、计算量较大 刚性方程、对稳定性要求高的场景、大步长计算
改进Euler方法 稳定 \((-1,0)\) 精度较高、显式无需迭代、计算效率高、兼顾精度与易用性 稳定性弱于梯形法 工程常规非刚性ODE求解、入门级数值计算、兼顾精度与效率的场景

表3 核心概念与结论总结

核心概念 定义 核心结论
局部截断误差 前一步计算值精确时，单步计算的误差 Euler方法：\(O(h^2)\)；梯形/改进Euler：\(O(h^3)\)
整体截断误差 递推到第\(m\)步时，数值解与精确解的总误差 Euler方法：\(O(h)\)；梯形/改进Euler：\(O(h^2)\)
收敛性 步长\(h \to 0\)时，数值解是否收敛于精确解 满足Lipschitz条件时，Euler系列方法均收敛
零稳定性 初始扰动是否会随计算步数无限放大 满足Lipschitz条件时，Euler系列方法均零稳定

Runge-Kutta（龙格-库塔）方法 完整讲解与推导证明

承接前文的Euler方法，我们已经知道单步法的核心是通过离散化将微分方程转化为递推格式。但Euler类方法精度有限，若要提升精度，直接使用Taylor展开法需要计算\(f(x,y)\)的高阶偏导数，计算量极大且编程实现复杂。Runge-Kutta（简称RK）方法正是为解决这一问题提出的：它通过多个节点处\(f(x,y)\)函数值的线性组合构造递推格式，完全避免了高阶求导，同时能达到任意指定的精度阶数，是目前工程中求解非刚性常微分方程最主流的单步数值方法。

一、显式Runge-Kutta方法的一般形式与阶条件

1.1 s级显式RK方法的核心定义

针对一阶常微分方程初值问题

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

s级显式Runge-Kutta方法的一般形式为：

\[y_{m+1} = y_m + h \sum_{i=1}^s b_i K_i \tag{1.3.2} \]

其中各阶段导数\(K_i\)的计算式为：

\[\begin{cases} K_1 = f(x_m, y_m) \\ K_i = f\left(x_m + c_i h,\ y_m + h \sum_{j=1}^{i-1} a_{ij} K_j\right), \quad i=2,3,\dots,s \end{cases} \tag{1.3.3} \]

参数含义

  • \(s\)：方法的级数，即每一步需要计算的\(f(x,y)\)的次数；
  • \(c_i\)：节点横坐标的偏移量，决定第\(i\)个导数\(K_i\)的计算位置；
  • \(a_{ij}\)：阶段导数的权重系数，决定\(K_j\)对\(K_i\)的贡献；
  • \(b_i\)：最终的组合权重，决定每个\(K_i\)对\(y_{m+1}\)的贡献；
  • 显式的核心特征：当\(j\geq i\)时，\(a_{ij}=0\)，即\(K_i\)仅依赖于前面已计算出的\(K_1,K_2,\dots,K_{i-1}\)，可依次递推计算，无需解方程组。

1.2 Butcher表：RK方法的标准表示

为了简洁清晰地表示RK方法的所有参数，我们使用Butcher表（布彻表），格式如下：

\[\begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \dots & a_{1s} \\ c_2 & a_{21} & a_{22} & \dots & a_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \dots & a_{ss} \\ \hline & b_1 & b_2 & \dots & b_s \end{array} \]

对于显式RK方法，\(a_{ij}=0\ (j\geq i)\)，因此Butcher表的\(a\)矩阵是严格下三角矩阵，例如显式Euler方法（一级一阶）的Butcher表为：

\[\begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \]

1.3 阶条件的完整推导（以二级二阶RK为例）

RK方法的核心是通过参数匹配，让方法的展开式与精确解的Taylor展开式尽可能多的项重合，从而达到指定的精度阶数。

阶数的定义

若方法的局部截断误差为\(O(h^{p+1})\)（即前一步计算值精确时，单步计算的误差为\(h^{p+1}\)阶），则称该方法是p阶的。要构造p阶RK方法，需要让方法的展开式与\(y(x_{m+1})\)的Taylor展开式的前\(p+1\)项完全相等，由此得到的参数方程组称为阶条件。

二级二阶RK方法的阶条件推导

我们以\(s=2\)（二级）的显式RK方法为例，完整推导二阶方法的阶条件。

二级显式RK的一般形式为：

\[\begin{cases} y_{m+1} = y_m + h\left(b_1 K_1 + b_2 K_2\right) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + c_2 h,\ y_m + h a_{21} K_1\right) \end{cases} \]

步骤1：精确解的Taylor展开

假设\(y_m = y(x_m)\)（局部截断误差的定义，前一步计算值精确），将\(y(x_{m+1})=y(x_m+h)\)在\(x_m\)处做Taylor展开到\(h^2\)项：

\[y(x_{m+1}) = y(x_m) + h y'(x_m) + \frac{h^2}{2} y''(x_m) + O(h^3) \]

由微分方程\(y'(x)=f(x,y(x))\)，得：

  • \(y'(x_m) = f(x_m, y_m) \triangleq f\)
  • \(y''(x_m) = \frac{d}{dx}f(x,y(x))\bigg|_{x=x_m} = f_x(x_m,y_m) + f_y(x_m,y_m) \cdot y'(x_m) = f_x + f_y f\)

因此精确解的展开式为：

\[y(x_{m+1}) = y_m + h f + \frac{h^2}{2}\left(f_x + f_y f\right) + O(h^3) \tag{A} \]

步骤2：RK方法的Taylor展开

我们将\(y_{m+1}\)也展开到\(h^2\)项。首先对\(K_2\)做二元Taylor展开：

\[K_2 = f\left(x_m + c_2 h,\ y_m + h a_{21} f\right) \]

二元函数的一阶Taylor展开为：

\[f(x+\Delta x, y+\Delta y) = f(x,y) + \Delta x f_x + \Delta y f_y + O(\Delta x^2 + \Delta y^2) \]

代入\(\Delta x = c_2 h\)，\(\Delta y = h a_{21} f\)，得：

\[K_2 = f + c_2 h f_x + h a_{21} f \cdot f_y + O(h^2) \]

将\(K_1,K_2\)代入\(y_{m+1}\)的表达式，展开整理：

\[\begin{align*} y_{m+1} &= y_m + h\left[b_1 f + b_2 \left(f + c_2 h f_x + h a_{21} f f_y + O(h^2)\right)\right] \\ &= y_m + h (b_1 + b_2) f + h^2 \left(b_2 c_2 f_x + b_2 a_{21} f f_y\right) + O(h^3) \tag{B} \end{align*} \]

步骤3：系数匹配得到阶条件

要让方法达到二阶精度，必须让(A)和(B)的\(h^0,h^1,h^2\)项系数完全相等，由此得到方程组：

  1. \(h^1\)项系数相等：\(\boldsymbol{b_1 + b_2 = 1}\)
  2. \(h^2\)项\(f_x\)系数相等：\(\boldsymbol{b_2 c_2 = \frac{1}{2}}\)
  3. \(h^2\)项\(f f_y\)系数相等：\(\boldsymbol{b_2 a_{21} = \frac{1}{2}}\)

这就是二级二阶RK方法的阶条件。3个方程包含4个未知参数，因此存在1个自由参数，取不同的自由参数即可得到不同的二级二阶RK格式。

常见的二级二阶RK格式

  1. 中点公式：取自由参数\(c_2=\frac{1}{2}\)，解得\(b_2=1,\ b_1=0,\ a_{21}=\frac{1}{2}\)，格式为：

    \[\begin{cases} y_{m+1} = y_m + h K_2 \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_1}{2}\right) \end{cases} \]

    Butcher表：

    \[\begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array} \]

  2. 改进Euler公式（Heun公式）：取自由参数\(c_2=1\)，解得\(b_2=\frac{1}{2},\ b_1=\frac{1}{2},\ a_{21}=1\)，格式为：

    \[\begin{cases} y_{m+1} = y_m + \frac{h}{2}(K_1 + K_2) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + h,\ y_m + h K_1\right) \end{cases} \]

    Butcher表：

    \[\begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array} \]

  3. Ralston公式：取自由参数\(c_2=\frac{2}{3}\)，解得\(b_2=\frac{3}{4},\ b_1=\frac{1}{4},\ a_{21}=\frac{2}{3}\)，格式为：

    \[\begin{cases} y_{m+1} = y_m + h\left(\frac{1}{4}K_1 + \frac{3}{4}K_2\right) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{2h}{3},\ y_m + \frac{2h K_1}{3}\right) \end{cases} \]

1.4 高阶RK方法的阶条件说明

随着精度阶数\(p\)的提升，阶条件的数量会快速增长，且存在Butcher壁垒：

  • 显式RK方法中，\(s\)级方法最多能达到\(s\)阶精度；
  • 阶条件数量：\(p=1\)需1个，\(p=2\)需2个，\(p=3\)需4个，\(p=4\)需8个，\(p=5\)需17个；
  • 当\(p\geq5\)时，达到\(p\)阶精度需要的级数\(s>p\)（例如五阶方法至少需要6级），计算量会显著提升。

因此工程中最常用的是四级四阶经典RK方法，完美平衡了精度与计算量。

二、常用的显式Runge-Kutta格式

2.1 三级三阶RK格式

（1）三级三阶Heun公式

\[\begin{cases} y_{m+1} = y_m + \frac{h}{4}\left(K_1 + 3K_3\right) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{h}{3},\ y_m + \frac{h K_1}{3}\right) \\ K_3 = f\left(x_m + \frac{2h}{3},\ y_m + \frac{2h K_2}{3}\right) \end{cases} \]

Butcher表：

\[\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array} \]

（2）三级三阶Kutta公式

\[\begin{cases} y_{m+1} = y_m + \frac{h}{6}\left(K_1 + 4K_2 + K_3\right) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_1}{2}\right) \\ K_3 = f\left(x_m + h,\ y_m - h K_1 + 2h K_2\right) \end{cases} \]

2.2 四级四阶经典RK公式（RK4）

这是工程中最常用的RK方法，局部截断误差为\(O(h^5)\)，整体截断误差为\(O(h^4)\)，格式如下：

\[\boldsymbol{ \begin{cases} y_{m+1} = y_m + \frac{h}{6}\left(K_1 + 2K_2 + 2K_3 + K_4\right) \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_1}{2}\right) \\ K_3 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_2}{2}\right) \\ K_4 = f\left(x_m + h,\ y_m + h K_3\right) \end{cases} } \]

Butcher表：

\[\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 2/6 & 2/6 & 1/6 \end{array} \]

2.3 其他常用高阶显式RK格式

  1. 四级四阶Gill公式：具有减小舍入误差的优势，适合长区间计算

    \[\begin{cases} y_{m+1} = y_m + \frac{h}{6}\left[K_1 + (2-\sqrt{2})K_2 + (2+\sqrt{2})K_3 + K_4\right] \\ K_1 = f(x_m, y_m) \\ K_2 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_1}{2}\right) \\ K_3 = f\left(x_m + \frac{h}{2},\ y_m + \frac{(\sqrt{2}-1)h K_1}{2} + \frac{(2-\sqrt{2})h K_2}{2}\right) \\ K_4 = f\left(x_m + h,\ y_m - \frac{\sqrt{2} h K_2}{2} + \frac{(2+\sqrt{2})h K_3}{2}\right) \end{cases} \]

  2. 六级五阶Nystrom公式、七级六阶Lawson公式：用于高精度计算，Butcher表见教材原文。

三、隐式Runge-Kutta方法

显式RK方法的绝对稳定区间有限，不适合求解刚性微分方程，且级数\(s\)最多达到\(s\)阶精度。隐式Runge-Kutta方法突破了这些限制，是求解刚性方程的核心方法。

3.1 隐式RK的定义与分类

定义

s级RK方法的一般形式为：

\[\begin{cases} y_{m+1} = y_m + h \sum_{i=1}^s b_i K_i \\ K_i = f\left(x_m + c_i h,\ y_m + h \sum_{j=1}^s a_{ij} K_j\right), \quad i=1,2,\dots,s \end{cases} \tag{1.3.11} \]

若存在\(i\leq j\)使得\(a_{ij}\neq0\)，则\(K_i\)的表达式中包含自身或后续的\(K_j\)，无法依次递推，需要求解非线性方程组，这类方法称为隐式Runge-Kutta方法。

分类

  1. 全隐式RK方法：\(a\)矩阵为满矩阵，精度最高，稳定性最好，每一步需要求解\(s\)维非线性方程组；
  2. 对角隐式RK（DIRK）：\(a_{ij}=0\ (j>i)\)，即\(a\)矩阵为下三角矩阵，对角元非零。每个\(K_i\)仅依赖于前面的\(K_1,\dots,K_i\)，每一步只需依次求解\(s\)个一维非线性方程，计算量远小于全隐式；
  3. 单对角隐式RK（SDIRK）：DIRK的对角元全部相等，进一步简化计算，是工程中最常用的隐式RK类型。

3.2 隐式RK的核心优势

  1. 精度上限更高：s级全隐式RK方法最高可达到\(2s\)阶精度（显式RK最多\(s\)阶），例如1级隐式中点公式为2阶，2级Gauss型RK为4阶；
  2. 稳定性极强：多数隐式RK方法是A-稳定的（绝对稳定区域包含整个复平面左半平面），甚至L-稳定，完美适配刚性微分方程的求解；
  3. 保结构特性：Gauss型隐式RK是辛算法，适合求解Hamilton系统等保结构问题。

3.3 隐式RK的收敛性证明

隐式RK方法的核心问题是：非线性方程组是否存在唯一解？我们通过压缩映射原理给出收敛性定理的完整证明。

定理1.3.1

设\(f(x,y)\)关于\(y\)满足Lipschitz条件，Lipschitz常数为\(L\)，若步长\(h\)满足：

\[h L \max_{1\leq i\leq s} \sum_{j=1}^s |a_{ij}| < 1 \tag{1.3.12} \]

则隐式RK方法的阶段导数\(K_1,K_2,\dots,K_s\)存在唯一解。

完整证明

将阶段导数写成向量形式\(\boldsymbol{K} = (K_1,K_2,\dots,K_s)^T\)，定义映射\(\boldsymbol{F}(\boldsymbol{K})\)，其中第\(i\)个分量为：

\[F_i(\boldsymbol{K}) = f\left(x_m + c_i h,\ y_m + h \sum_{j=1}^s a_{ij} K_j\right), \quad i=1,\dots,s \]

则隐式RK的阶段导数满足不动点方程\(\boldsymbol{K} = \boldsymbol{F}(\boldsymbol{K})\)。我们只需证明\(\boldsymbol{F}\)是压缩映射，即可由压缩映射原理得到不动点的存在唯一性。

对任意两个向量\(\boldsymbol{K}\)和\(\boldsymbol{K}'\)，计算分量的差值：

\[\begin{align*} |F_i(\boldsymbol{K}) - F_i(\boldsymbol{K}')| &= \left| f\left(x_m + c_i h,\ y_m + h \sum_{j=1}^s a_{ij} K_j\right) - f\left(x_m + c_i h,\ y_m + h \sum_{j=1}^s a_{ij} K'_j\right) \right| \\ &\leq L \cdot \left| h \sum_{j=1}^s a_{ij} (K_j - K'_j) \right| \\ &\leq h L \sum_{j=1}^s |a_{ij}| \cdot |K_j - K'_j| \end{align*} \]

取无穷范数\(\|\boldsymbol{K}\|_\infty = \max_{1\leq i\leq s} |K_i|\)，则：

\[\|\boldsymbol{F}(\boldsymbol{K}) - \boldsymbol{F}(\boldsymbol{K}')\|_\infty \leq h L \left( \max_{1\leq i\leq s} \sum_{j=1}^s |a_{ij}| \right) \cdot \|\boldsymbol{K} - \boldsymbol{K}'\|_\infty \]

当\(h L \max_{i} \sum_j |a_{ij}| < 1\)时，\(\boldsymbol{F}\)是压缩映射，根据压缩映射原理，存在唯一的不动点\(\boldsymbol{K}\)，即阶段导数存在唯一解。定理得证。

3.4 常用隐式RK格式

  1. 一级二阶隐式中点公式（最基础的隐式RK，A-稳定）

    \[\begin{cases} y_{m+1} = y_m + h K_1 \\ K_1 = f\left(x_m + \frac{h}{2},\ y_m + \frac{h K_1}{2}\right) \end{cases} \]

    Butcher表：

    \[\begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array} \]

  2. 二级四阶Gauss型RK公式（A-稳定、L-稳定，辛算法）
    Butcher表：

    \[\begin{array}{c|cc} \frac{3-\sqrt{3}}{6} & \frac{1}{4} & \frac{3-2\sqrt{3}}{12} \\ \frac{3+\sqrt{3}}{6} & \frac{3+2\sqrt{3}}{12} & \frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \]

  3. 二级二阶Radau IIA公式、三级四阶Lobatto IIIA公式：用于边值问题和刚性方程求解，Butcher表见教材原文。

四、单步法的收敛性与稳定性理论

收敛性与稳定性是数值方法的核心理论基础，决定了方法是否可靠、是否能得到有意义的数值解。

4.1 单步法的一般形式与基本定义

所有单步法（Euler方法、RK方法）都可以写成统一的形式：

\[y_{m+1} = y_m + h \Phi(x_m, y_m, h) \tag{1.3.19} \]

其中\(\Phi(x,y,h)\)称为增量函数，例如：

  • 显式Euler方法：\(\Phi(x,y,h) = f(x,y)\)
  • 经典四阶RK方法：\(\Phi(x,y,h) = \frac{1}{6}(K_1 + 2K_2 + 2K_3 + K_4)\)

核心定义

  1. 局部截断误差：假设前一步计算值精确（\(y_m = y(x_m)\)），单步计算的误差

    \[T_m = y(x_{m+1}) - \left[ y(x_m) + h \Phi(x_m, y(x_m), h) \right] \]

  2. 相容阶：若\(T_m = O(h^{p+1})\)，则称方法是p阶相容的；当\(p\geq1\)时，称方法与微分方程相容，相容的充要条件是\(\Phi(x,y,0) = f(x,y)\)。
  3. 零稳定性：存在常数\(C>0\)和\(h_0>0\)，使得对任意\(0<h<h_0\)、任意两个数值解\(y_m\)和\(z_m\)，有

    \[|y_m - z_m| \leq C |y_0 - z_0| \]

    即初始值的微小扰动不会随计算步数无限放大。
  4. 收敛性：当步长\(h\to0\)时，\(\max_{0\leq m\leq N} |y(x_m) - y_m| \to 0\)，称方法是收敛的。

4.2 Lax等价定理：收敛性的充要条件

定理内容

对于单步法，方法收敛的充要条件是方法相容且零稳定。

完整证明

我们仅证明核心的充分性（相容+稳定⇒收敛），必要性可由收敛的定义直接推导。

充分性证明：
设整体截断误差\(e_m = y(x_m) - y_m\)，由局部截断误差的定义：

\[y(x_{m+1}) = y(x_m) + h \Phi(x_m, y(x_m), h) + T_m \]

数值解满足：

\[y_{m+1} = y_m + h \Phi(x_m, y_m, h) \]

两式相减得误差递推式：

\[e_{m+1} = e_m + h \left[ \Phi(x_m, y(x_m), h) - \Phi(x_m, y_m, h) \right] + T_m \]

两边取绝对值，假设增量函数\(\Phi\)关于\(y\)满足Lipschitz条件（RK方法、Euler方法均满足，因\(f\)满足Lipschitz条件），即存在\(L_\Phi>0\)，使得：

\[|\Phi(x,y,h) - \Phi(x,z,h)| \leq L_\Phi |y-z| \]

代入得：

\[|e_{m+1}| \leq (1 + h L_\Phi) |e_m| + |T_m| \]

方法是p阶相容的，因此存在常数\(C_T>0\)，使得\(|T_m| \leq C_T h^{p+1}\)，代入递推式：

\[|e_{m+1}| \leq (1 + h L_\Phi) |e_m| + C_T h^{p+1} \]

对该不等式递推展开：

\[|e_m| \leq (1 + h L_\Phi)^m |e_0| + C_T h^{p+1} \sum_{k=0}^{m-1} (1 + h L_\Phi)^k \]

等比数列求和得：

\[\sum_{k=0}^{m-1} (1 + h L_\Phi)^k = \frac{(1 + h L_\Phi)^m - 1}{h L_\Phi} \]

利用不等式\(1 + t \leq e^t\)（\(t>0\)），得\((1 + h L_\Phi)^m \leq e^{m h L_\Phi} \leq e^{L_\Phi (b-a)}\)（因\(mh \leq b-a\)），代入得：

\[|e_m| \leq e^{L_\Phi (b-a)} |e_0| + \frac{C_T}{L_\Phi} \left( e^{L_\Phi (b-a)} - 1 \right) h^p \]

若初始值精确（\(e_0=0\)），则\(|e_m| = O(h^p)\)，当\(h\to0\)时，\(|e_m|\to0\)，方法收敛。充分性得证。

4.3 绝对稳定性与绝对稳定区间

零稳定性是\(h\to0\)时的渐进稳定性，而实际计算中\(h\)是有限值，我们需要分析有限步长下方法的稳定性，即绝对稳定性。

模型方程与稳定性函数

采用线性测试方程（模型方程）：

\[y' = \lambda y, \quad \text{Re}(\lambda) < 0 \]

其中\(\lambda\)为复常数，Re(\(\lambda\))<0保证方程的精确解\(y(x)=y_0 e^{\lambda x}\)是衰减的。

将单步法代入模型方程，得到递推式：

\[y_{m+1} = R(z) y_m, \quad z = \lambda h \]

其中\(R(z)\)称为方法的稳定性函数。

绝对稳定的定义

若\(|R(z)| < 1\)，则称方法在\(z\)点是绝对稳定的；复平面上所有满足\(|R(z)|<1\)的\(z\)的集合称为绝对稳定区域，绝对稳定区域与实轴的交集称为绝对稳定区间。

常用RK方法的稳定性函数与绝对稳定区间

方法 稳定性函数\(R(z)\) 实轴绝对稳定区间
显式Euler（一级一阶） \(1+z\) \((-2, 0)\)
二级二阶中点公式 \(1+z+\frac{z^2}{2}\) \((-2, 0)\)
三级三阶Kutta公式 \(1+z+\frac{z^2}{2}+\frac{z^3}{6}\) \((-2.51, 0)\)
四级四阶经典RK \(1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}\) \((-2.78, 0)\)
隐式Euler（一级一阶） \(\frac{1}{1-z}\) \((-\infty, 0)\)（A-稳定）
隐式中点公式（一级二阶） \(\frac{1+z/2}{1-z/2}\) \((-\infty, 0)\)（A-稳定）

核心结论：显式RK方法的绝对稳定区间都是有限的，步长不能太大；隐式RK方法可实现A-稳定，步长不受稳定性限制，仅受精度限制，完美适配刚性方程。

五、核心知识点系统性归纳总结

表1 常用显式Runge-Kutta格式对比

方法名称 级数 阶数 每步\(f\)计算次数 局部截断误差 绝对稳定区间 核心特点
显式Euler方法 1 1 1次 \(O(h^2)\) \((-2,0)\) 结构最简单、计算量最小、精度最低
二级二阶中点公式 2 2 2次 \(O(h^3)\) \((-2,0)\) 单步精度高、相位误差小
改进Euler公式 2 2 2次 \(O(h^3)\) \((-2,0)\) 直观易懂、工程入门常用
三级三阶Heun公式 3 3 3次 \(O(h^4)\) \((-2.51,0)\) 平衡精度与计算量
四级四阶经典RK公式 4 4 4次 \(O(h^5)\) \((-2.78,0)\) 工程最常用、精度与计算量完美平衡
四级四阶Gill公式 4 4 4次 \(O(h^5)\) \((-2.78,0)\) 舍入误差小、适合长区间计算

表2 隐式Runge-Kutta格式对比

方法名称 级数 阶数 稳定性 核心特点 适用场景
隐式中点公式 1 2 A-稳定、L-稳定 结构最简单、辛算法 刚性方程、Hamilton系统
二级四阶Gauss型RK 2 4 A-稳定、L-稳定 高精度、辛算法 高精度刚性方程、保结构计算
对角隐式RK（DIRK） s s A-稳定 计算量小于全隐式、易实现 工程刚性方程求解
Radau IIA公式 s 2s-1 A-稳定、L-稳定 适合边值问题 微分方程边值问题、刚性方程

表3 单步法核心理论结论总结

核心概念 核心定义 关键结论
相容性 方法的增量函数满足\(\Phi(x,y,0)=f(x,y)\) 方法收敛的必要条件，对应局部截断误差\(O(h^2)\)
零稳定性 初始扰动不会随步数无限放大 方法收敛的必要条件，所有RK方法均满足零稳定性
收敛性 \(h\to0\)时数值解收敛于精确解 充要条件：相容+零稳定（Lax等价定理）
绝对稳定性 有限步长下误差不会指数增长 显式RK绝对稳定区间有限，隐式RK可实现A-稳定
Butcher壁垒 显式RK的级数与阶数的关系 s级显式RK最多达到s阶精度，s级隐式RK最高可达2s阶

线性多步法 完整讲解与推导证明

承接前文的单步法（Euler、Runge-Kutta方法），我们已经知道单步法的核心是仅利用前一个节点的信息计算下一个节点值，若要提升精度，需在每一步多次计算\(f(x,y)\)的值，计算成本随阶数快速上升。线性多步法突破了这一限制：它充分利用之前多个节点已经计算出的\(y\)值与\(f\)值，构造递推格式，每一步仅需计算1次新的\(f(x,y)\)值，即可达到高阶精度，是长区间、大规模常微分方程求解中效率极高的数值方法。

一、线性多步法的基本定义与核心概念

1.1 k步线性多步法的一般形式

针对一阶常微分方程初值问题

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

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

\[\boldsymbol{y_{m+k} + \sum_{j=0}^{k-1} \alpha_j y_{m+j} = h \sum_{j=0}^k \beta_j f_{m+j}} \tag{1.4.1} \]

其中：

  • \(f_{m+j} = f(x_{m+j}, y_{m+j})\)，为节点\(x_{m+j}\)处的导数值；
  • \(x_{m+j} = x_0 + (m+j)h\)，等距步长\(h\)；
  • \(\alpha_j, \beta_j\)为与\(m\)无关的常数，且\(\alpha_0\)与\(\beta_0\)不同时为0；
  • 步数k：计算\(y_{m+k}\)时，需要用到前k个节点的\(y\)值\(y_m, y_{m+1}, \dots, y_{m+k-1}\)，因此称为k步法；
  • 显式与隐式的区分：
    • 若\(\beta_k = 0\)，等式右端不含\(f_{m+k}\)（即不含未知量\(y_{m+k}\)），可直接计算\(y_{m+k}\)，称为显式线性多步法；
    • 若\(\beta_k \neq 0\)，等式右端含\(f_{m+k}\)，未知量\(y_{m+k}\)出现在等式两端，需迭代求解，称为隐式线性多步法。

关键特性：k步法无法自启动，必须先用单步法（如Runge-Kutta方法）计算出前k个初始值\(y_0,y_1,\dots,y_{k-1}\)，才能启动递推计算。

1.2 局部截断误差、阶数与相容性

1.2.1 局部截断误差的定义

线性多步法的局部截断误差，是衡量单步计算精度的核心指标，定义为：
假设前k步的计算值均为精确值，即\(y_{m+j} = y(x_{m+j}),\ j=0,1,\dots,k-1\)，则用精确解代入方法格式得到的残差为局部截断误差\(R_{m,k}\)：

\[\boldsymbol{R_{m,k} = y(x_{m+k}) + \sum_{j=0}^{k-1} \alpha_j y(x_{m+j}) - h \sum_{j=0}^k \beta_j y'(x_{m+j})} \tag{1.4.2} \]

其中\(y'(x_{m+j}) = f(x_{m+j}, y(x_{m+j}))\)为精确解的导数。

1.2.2 阶数与相容性

  • p阶方法：若局部截断误差满足\(R_{m,k} = O(h^{p+1})\)，且\(p\geq1\)，则称该k步线性多步法是p阶精度的；
  • 相容性：若方法的阶数\(p\geq1\)，则称该方法与原微分方程相容，相容性是方法收敛的必要条件。

1.3 线性多步法的阶条件（Taylor展开推导）

我们通过Taylor展开，推导线性多步法达到p阶精度需要满足的条件（阶条件）。

将\(y(x_{m+j})\)和\(y'(x_{m+j})\)在\(x_m\)处做Taylor展开：

\[y(x_{m+j}) = y(x_m + jh) = \sum_{n=0}^\infty \frac{(jh)^n}{n!} y^{(n)}(x_m) \]

\[y'(x_{m+j}) = \sum_{n=0}^\infty \frac{(jh)^n}{n!} y^{(n+1)}(x_m) \]

将其代入局部截断误差的表达式(1.4.2)，整理得：

\[\begin{align*} R_{m,k} &= \sum_{n=0}^\infty \frac{(kh)^n}{n!} y^{(n)}(x_m) + \sum_{j=0}^{k-1} \alpha_j \sum_{n=0}^\infty \frac{(jh)^n}{n!} y^{(n)}(x_m) - h \sum_{j=0}^k \beta_j \sum_{n=0}^\infty \frac{(jh)^n}{n!} y^{(n+1)}(x_m) \\ &= \sum_{n=0}^\infty \frac{h^n}{n!} \left[ k^n + \sum_{j=0}^{k-1} \alpha_j j^n - n \sum_{j=0}^k \beta_j j^{n-1} \right] y^{(n)}(x_m) \end{align*} \]

（注：\(n=0\)时，\(n\sum\beta_j j^{n-1}=0\)）

要使方法达到p阶精度，即\(R_{m,k}=O(h^{p+1})\)，必须让\(h^0,h^1,\dots,h^p\)项的系数全部为0，由此得到p阶线性多步法的阶条件方程组：

\[\boldsymbol{ \begin{cases} C_0 = 1 + \sum_{j=0}^{k-1} \alpha_j = 0 \\ C_n = \frac{1}{n!} \left( k^n + \sum_{j=0}^{k-1} \alpha_j j^n - n \sum_{j=0}^k \beta_j j^{n-1} \right) = 0, \quad n=1,2,\dots,p \end{cases} } \]

其中\(C_{p+1}\)称为误差常数，局部截断误差可表示为\(R_{m,k} = C_{p+1} h^{p+1} y^{(p+1)}(x_m) + O(h^{p+2})\)。

核心结论：k步线性多步法最多可达到2k阶精度（隐式），显式k步线性多步法最多可达到k阶精度。

二、Adams线性多步法（最常用的线性多步法）

Adams方法是工程中应用最广泛的线性多步法，其核心思想是从微分方程的积分形式出发，通过插值多项式近似被积函数，构造递推格式，分为显式Adams-Bashforth方法（外插法）和隐式Adams-Moulton方法（内插法）两大类。

2.1 Adams方法的核心出发点

对微分方程\(y'(s)=f(s,y(s))\)在区间\([x_m, x_{m+1}]\)上积分，得到精确等式：

\[\boldsymbol{y(x_{m+1}) = y(x_m) + \int_{x_m}^{x_{m+1}} f(s,y(s)) ds} \tag{1.4.4} \]

Adams方法的核心是：用节点处的\(f\)值构造插值多项式\(P(s)\)，近似代替被积函数\(f(s,y(s))\)，通过积分\(\int_{x_m}^{x_{m+1}} P(s) ds\)得到递推格式。

2.2 显式Adams-Bashforth方法（外插法）

2.2.1 格式推导

显式Adams方法采用牛顿向后插值多项式，插值节点取积分区间左侧的k个节点：\(x_m, x_{m-1}, \dots, x_{m-k+1}\)，共k个节点，构造k-1次插值多项式\(P_{k-1}(s)\)，近似\(f(s,y(s))\)。

引入向后差分算子\(\nabla\)：\(\nabla f_m = f_m - f_{m-1}\)，\(\nabla^n f_m = \nabla^{n-1} f_m - \nabla^{n-1} f_{m-1}\)，牛顿向后插值公式为：

\[f(x_m + \tau h) = \sum_{j=0}^{k-1} (-1)^j \binom{-\tau}{j} \nabla^j f_m + R_{k-1}(x_m+\tau h) \]

其中\(\binom{-\tau}{j} = \frac{\tau(\tau+1)\dots(\tau+j-1)}{j!}\)，余项\(R_{k-1}(s) = h^k \binom{-\tau}{k} y^{(k+1)}(\xi),\ \xi\in(x_{m-k+1},x_m)\)。

将插值多项式代入积分等式(1.4.4)，令\(s = x_m + \tau h\)，则\(ds = h d\tau\)，积分区间\(\tau\in[0,1]\)，得：

\[y(x_{m+1}) \approx y(x_m) + h \int_0^1 \sum_{j=0}^{k-1} (-1)^j \binom{-\tau}{j} \nabla^j f_m d\tau \]

交换求和与积分顺序，定义系数\(a_j = \int_0^1 (-1)^j \binom{-\tau}{j} d\tau\)，则得到k步显式Adams-Bashforth方法的一般形式：

\[\boldsymbol{y_{m+1} = y_m + h \sum_{j=0}^{k-1} a_j \nabla^j f_m} \tag{1.4.10} \]

将向后差分展开为\(f\)的显式形式，可得到更常用的格式：

\[\boldsymbol{y_{m+1} = y_m + h \sum_{j=0}^{k-1} b_{k,j} f_{m-j}} \]

其中\(b_{k,j}\)为显式系数，可通过\(a_j\)推导得到。

2.2.2 系数与局部截断误差

  • 阶数：k步显式Adams-Bashforth方法为k阶精度，局部截断误差为\(R_{m,k} = C_{k+1} h^{k+1} y^{(k+1)}(x_m) + O(h^{k+2})\)；
  • 系数\(a_j\)的递推公式：通过生成函数可推导出\(a_j\)满足\(a_0=1\)，\(a_j + \frac{1}{2}a_{j-1} + \frac{1}{3}a_{j-2} + \dots + \frac{1}{j+1}a_0 = 1\)，\(j\geq1\)；
  • 常用的\(a_j\)数值表：
j 0 1 2 3 4 5
\(a_j\) 1 1/2 5/12 3/8 251/720 95/288

2.2.3 常用的显式Adams-Bashforth格式

  1. 1步1阶（显式Euler方法）：

    \[y_{m+1} = y_m + h f_m \]

    局部截断误差：\(R_m = \frac{1}{2}h^2 y''(x_m) + O(h^3)\)

  2. 2步2阶：

    \[y_{m+1} = y_m + \frac{h}{2}\left(3f_m - f_{m-1}\right) \]

    局部截断误差：\(R_m = \frac{5}{12}h^3 y'''(x_m) + O(h^4)\)

  3. 3步3阶：

    \[y_{m+1} = y_m + \frac{h}{12}\left(23f_m - 16f_{m-1} + 5f_{m-2}\right) \]

    局部截断误差：\(R_m = \frac{3}{8}h^4 y^{(4)}(x_m) + O(h^5)\)

  4. 4步4阶（最4步4阶（最常用）：

    \[\boldsymbol{y_{m+1} = y_m + \frac{h}{24}\left(55f_m - 59f_{m-1} + 37f_{m-2} - 9f_{m-3}\right)} \]

    局部截断误差：\(R_m = \frac{251}{720}h^5 y^{(5)}(x_m) + O(h^6)\)

2.3 隐式Adams-Moulton方法（内插法）

2.3.1 格式推导

隐式Adams方法的插值节点包含积分区间的右端点\(x_{m+1}\)，取\(x_{m+1},x_m,\dots,x_{m-k+2}\)共k个节点，构造k-1次牛顿向后插值多项式，插值区间包含积分区间\([x_m,x_{m+1}]\)，因此称为内插法，得到的是隐式格式。

同样代入积分等式(1.4.4)，积分后得到k步隐式Adams-Moulton方法的一般形式：

\[\boldsymbol{y_{m+1} = y_m + h \sum_{j=0}^{k-1} a_j^* \nabla^j f_{m+1}} \]

展开为显式\(f\)的形式：

\[\boldsymbol{y_{m+1} = y_m + h \sum_{j=0}^{k-1} b_{k,j}^* f_{m+1-j}} \]

其中\(b_{k,0}^* \neq 0\)，因此等式右端含\(f_{m+1}=f(x_{m+1},y_{m+1})\)，为隐式格式。

2.3.2 系数与局部截断误差

  • 阶数：k步隐式Adams-Moulton方法为k阶精度，比同步数的显式Adams方法精度高一阶，局部截断误差为\(O(h^{k+1})\)；
  • 系数\(a_j^*\)满足递推：\(a_0^*=1\)，\(a_j^* + \frac{1}{2}a_{j-1}^* + \dots + \frac{1}{j+1}a_0^* = 0\)，\(j\geq1\)；
  • 常用的\(a_j^*\)数值表：
j 0 1 2 3 4
\(a_j^*\) 1 -1/2 -1/12 -1/24 -19/720

2.3.3 常用的隐式Adams-Moulton格式

  1. 1步2阶（梯形方法）：

    \[y_{m+1} = y_m + \frac{h}{2}\left(f_{m+1} + f_m\right) \]

    局部截断误差：\(R_m = -\frac{1}{12}h^3 y'''(x_m) + O(h^4)\)

  2. 2步3阶：

    \[y_{m+1} = y_m + \frac{h}{12}\left(5f_{m+1} + 8f_m - f_{m-1}\right) \]

    局部截断误差：\(R_m = -\frac{1}{24}h^4 y^{(4)}(x_m) + O(h^5)\)

  3. 3步4阶（最常用）：

    \[\boldsymbol{y_{m+1} = y_m + \frac{h}{24}\left(9f_{m+1} + 19f_m - 5f_{m-1} + f_{m-2}\right)} \]

    局部截断误差：\(R_m = -\frac{19}{720}h^5 y^{(5)}(x_m) + O(h^6)\)

  4. 4步5阶：

    \[y_{m+1} = y_m + \frac{h}{720}\left(251f_{m+1} + 646f_m - 264f_{m-1} + 106f_{m-2} - 19f_{m-3}\right) \]

2.3.4 隐式Adams方法的核心优势

  1. 精度更高：同步数下，隐式Adams比显式Adams精度高一阶；
  2. 误差常数更小：同阶数下，隐式方法的误差常数远小于显式方法，单步精度更高；
  3. 稳定性更好：隐式Adams的绝对稳定区间远大于显式Adams，甚至部分格式是A-稳定的，适合求解刚性方程。

2.4 Adams预测-校正格式

为了兼顾显式方法的计算效率和隐式方法的精度与稳定性，工程中最常用的是Adams预测-校正格式，核心思路是：

  1. 预测步：用显式Adams-Bashforth方法计算\(y_{m+1}\)的预测值\(\tilde{y}_{m+1}\)；
  2. 校正步：将预测值代入隐式Adams-Moulton方法的右端，计算校正值\(y_{m+1}\)，仅做一次校正，无需迭代。

最常用的4阶Adams预测-校正格式

\[\begin{cases} \text{预测步（4步4阶显式Adams）：} \quad \tilde{y}_{m+1} = y_m + \frac{h}{24}\left(55f_m - 59f_{m-1} + 37f_{m-2} - 9f_{m-3}\right) \\ \text{计算预测导数值：} \quad \tilde{f}_{m+1} = f(x_{m+1}, \tilde{y}_{m+1}) \\ \text{校正步（3步4阶隐式Adams）：} \quad y_{m+1} = y_m + \frac{h}{24}\left(9\tilde{f}_{m+1} + 19f_m - 5f_{m-1} + f_{m-2}\right) \\ \text{更新导数值：} \quad f_{m+1} = f(x_{m+1}, y_{m+1}) \end{cases} \]

核心特性：该格式为4阶精度，每步仅需计算1次新的\(f\)值（预测步的\(f\)均为历史值，校正步仅需计算\(\tilde{f}_{m+1}\)），计算效率远高于同阶的Runge-Kutta方法，是长区间非刚性ODE求解的首选方法。

三、其他常用线性多步法

除了Adams方法，工程中还有几类常用的线性多步法，我们通过阶条件推导其格式与特性。

3.1 Milne方法（4阶显式多步法）

Milne方法是基于Simpson数值积分公式构造的4阶显式4步方法，格式为：

\[\boldsymbol{y_{m+4} = y_m + \frac{4h}{3}\left(2f_{m+3} - f_{m+2} + 2f_{m+1}\right)} \]

  • 局部截断误差：\(R_m = \frac{14}{45}h^5 y^{(5)}(x_m) + O(h^6)\)，4阶精度；
  • 特点：精度高，格式简单，但零稳定性较差，存在寄生根，容易出现数值振荡。

3.2 Simpson方法（4阶隐式多步法）

基于Simpson积分公式构造的2步4阶隐式方法，格式为：

\[\boldsymbol{y_{m+2} = y_m + \frac{h}{3}\left(f_{m+2} + 4f_{m+1} + f_m\right)} \]

  • 局部截断误差：\(R_m = -\frac{1}{90}h^5 y^{(5)}(x_m) + O(h^6)\)，4阶精度；
  • 特点：精度极高，误差常数极小，但零稳定性差，存在模为1的重根，仅适合求解边值问题，不适合初值问题的长区间计算。

3.3 Hamming方法（4阶隐式多步法）

Hamming方法是为了改进Milne方法的稳定性构造的4阶隐式4步方法，格式为：

\[\boldsymbol{y_{m+4} = \frac{1}{8}\left(9y_{m+3} - y_{m+1}\right) + \frac{3h}{8}\left(f_{m+4} + 2f_{m+3} - 2f_{m+1} - f_m\right)} \]

  • 局部截断误差：\(R_m = -\frac{1}{40}h^5 y^{(5)}(x_m) + O(h^6)\)，4阶精度；
  • 特点：零稳定性好，无寄生根的数值振荡，是工程中替代Milne方法的首选，常与Milne方法结合构成预测-校正格式。

四、线性多步法的收敛性与稳定性理论

4.1 零稳定性（根条件）

线性多步法的零稳定性，是衡量初始扰动是否会随计算步数无限放大的核心指标，也是方法收敛的必要条件。

4.1.1 第一特征多项式与根条件

对于k步线性多步法(1.4.1)，定义第一特征多项式：

\[\rho(\lambda) = \lambda^k + \sum_{j=0}^{k-1} \alpha_j \lambda^j \]

根条件：若\(\rho(\lambda)\)的所有根的模都不超过1，且模为1的根都是单根，则称该线性多步法满足根条件。

4.1.2 零稳定性的定义

称线性多步法是零稳定的，当且仅当该方法满足根条件。

核心结论：线性多步法的根条件是零稳定的充要条件。例如：

  • 显式Euler方法：\(\rho(\lambda)=\lambda-1\)，根为\(\lambda=1\)，满足根条件，零稳定；
  • Simpson方法：\(\rho(\lambda)=\lambda^2-1\)，根为\(\lambda=1\)和\(\lambda=-1\)，均为单根，模为1，满足根条件，但存在寄生根\(\lambda=-1\)，长区间计算会出现振荡；
  • 若方法存在模大于1的根，则不满足根条件，零不稳定。

4.2 收敛性：Lax等价定理

对于线性多步法，方法收敛的充要条件是：方法相容且零稳定。

  • 相容性：方法的阶数\(p\geq1\)，即满足\(C_0=0,C_1=0\)；
  • 零稳定性：满足根条件；
  • 收敛阶：若方法是p阶相容且零稳定的，则整体截断误差为\(O(h^p)\)，即p阶收敛。

4.3 绝对稳定性

与单步法类似，我们用线性测试方程\(y'=\lambda y\)（\(\text{Re}(\lambda)<0\)）分析线性多步法的绝对稳定性。

将测试方程代入k步线性多步法，得到线性差分方程：

\[y_{m+k} + \sum_{j=0}^{k-1} \alpha_j y_{m+j} = \lambda h \sum_{j=0}^k \beta_j y_{m+j} \]

令\(z=\lambda h\)，定义稳定性多项式：

\[\pi(\lambda,z) = \rho(\lambda) - z \sigma(\lambda) \]

其中\(\sigma(\lambda) = \sum_{j=0}^k \beta_j \lambda^j\)为第二特征多项式。

绝对稳定的定义

若对于给定的\(z\)，稳定性多项式\(\pi(\lambda,z)=0\)的所有根的模都小于1，则称方法在\(z\)点是绝对稳定的；复平面上所有满足该条件的\(z\)的集合称为绝对稳定区域。

常用线性多步法的绝对稳定区间（实轴）

方法 步数 阶数 实轴绝对稳定区间
显式Adams-Bashforth 1 1 (-2, 0)
显式Adams-Bashforth 2 2 (-1, 0)
显式Adams-Bashforth 3 3 (-0.55, 0)
显式Adams-Bashforth 4 4 (-0.30, 0)
隐式Adams-Moulton 1 2 (-∞, 0)
隐式Adams-Moulton 2 3 (-6, 0)
隐式Adams-Moulton 3 4 (-3, 0)
隐式Adams-Moulton 4 5 (-1.8, 0)

核心结论：同阶数下，隐式Adams方法的绝对稳定区间远大于显式Adams方法；显式Adams方法的绝对稳定区间随阶数升高而缩小，步长限制更严格。

五、核心知识点系统性归纳总结

表1 常用线性多步法格式与特性对比

方法名称 类型 步数 阶数 递推格式 局部截断误差主项 实轴绝对稳定区间 核心特点
显式Euler 显式 1 1 \(y_{m+1}=y_m+hf_m\) \(\frac{1}{2}h^2 y''(x_m)\) (-2,0) 结构最简单、可自启动、精度最低
2步2阶Adams-Bashforth 显式 2 2 \(y_{m+1}=y_m+\frac{h}{2}(3f_m-f_{m-1})\) \(\frac{5}{12}h^3 y'''(x_m)\) (-1,0) 2阶精度、显式、不可自启动
4步4阶Adams-Bashforth 显式 4 4 \(y_{m+1}=y_m+\frac{h}{24}(55f_m-59f_{m-1}+37f_{m-2}-9f_{m-3})\) \(\frac{251}{720}h^5 y^{(5)}(x_m)\) (-0.30,0) 最常用显式多步法、4阶精度、每步仅算1次f
梯形法 隐式 1 2 \(y_{m+1}=y_m+\frac{h}{2}(f_{m+1}+f_m)\) \(-\frac{1}{12}h^3 y'''(x_m)\) (-∞,0) A-稳定、2阶精度、隐式、可自启动
3步4阶Adams-Moulton 隐式 3 4 \(y_{m+1}=y_m+\frac{h}{24}(9f_{m+1}+19f_m-5f_{m-1}+f_{m-2})\) \(-\frac{19}{720}h^5 y^{(5)}(x_m)\) (-3,0) 最常用隐式多步法、4阶精度、稳定性好
4阶Adams预测-校正 显式+隐式 4 4 预测：4阶显式Adams，校正：4阶隐式Adams \(O(h^5)\) (-0.7,0) 兼顾效率与精度、每步仅算1次f、工程首选
Milne方法 显式 4 4 \(y_{m+4}=y_m+\frac{4h}{3}(2f_{m+3}-f_{m+2}+2f_{m+1})\) \(\frac{14}{45}h^5 y^{(5)}(x_m)\) 不稳定区间大 精度高、稳定性差、易振荡
Hamming方法 隐式 4 4 \(y_{m+4}=\frac{1}{8}(9y_{m+3}-y_{m+1})+\frac{3h}{8}(f_{m+4}+2f_{m+3}-2f_{m+1}-f_m)\) \(-\frac{1}{40}h^5 y^{(5)}(x_m)\) (-0.7,0) 稳定性好、无振荡、替代Milne方法

表2 线性多步法核心理论结论总结

核心概念 核心定义 关键结论
相容性 方法阶数\(p\geq1\)，满足\(C_0=0,C_1=0\) 方法收敛的必要条件
零稳定性 方法满足根条件（第一特征多项式的根模≤1，模1的根为单根） 方法收敛的必要条件，线性多步法的核心稳定性要求
收敛性 \(h\to0\)时数值解一致收敛于精确解 充要条件：相容+零稳定（Lax等价定理）
绝对稳定性 稳定性多项式的所有根模<1 隐式方法的绝对稳定区间远大于显式方法，步长限制更宽松
阶数上限 k步方法的最高可达精度 显式k步法最高k阶，隐式k步法最高2k阶

表3 线性多步法与单步法（Runge-Kutta）核心对比

特性 线性多步法 单步法（Runge-Kutta）
自启动 不可自启动，需单步法提供初始值 可自启动，仅需初始值\(y_0\)
计算效率 同阶精度下，每步仅需计算1次f，效率极高 同阶精度下，每步需多次计算f，效率较低
步长调整 步长调整复杂，需重新计算初始值 步长调整灵活，可随时改变步长
精度上限 隐式多步法可达到2k阶，上限更高 显式RK最多s阶，上限较低
稳定性 隐式多步法稳定性好，适合刚性方程 显式RK绝对稳定区间有限，不适合刚性方程
适用场景 长区间、大规模、非刚性/刚性方程求解 短区间、步长需频繁调整的场景、提供初始值

线性多步法的稳定性、收敛性与局部截断误差 完整讲解

本节是线性多步法的核心理论部分，系统阐述了线性多步法的零稳定性、收敛性、局部截断误差与绝对稳定性，是判断线性多步法是否可靠、能否用于实际计算的理论基石。

一、线性多步法的差分方程与特征方程

1.1 线性多步法的一般形式

k步线性多步法的标准形式为：

\[\boldsymbol{y_{m+k} + \sum_{j=0}^{k-1} \alpha_j y_{m+j} = h \sum_{j=0}^k \beta_j f_{m+j}} \tag{1.4.1} \]

其中 \(f_{m+j} = f(x_{m+j}, y_{m+j})\)，\(\alpha_j, \beta_j\) 为常数，且 \(\alpha_0, \beta_0\) 不同时为0。

1.2 对应的线性齐次差分方程

将方法中的 \(y_{m+j}\) 视为未知序列，令 \(h=0\)（或忽略右端项），得到对应的线性齐次差分方程：

\[\boldsymbol{y_{m+k} + \sum_{j=0}^{k-1} \alpha_j y_{m+j} = 0} \tag{1.5.1} \]

这是分析线性多步法零稳定性的核心方程。

1.3 特征方程与特征根

为求解差分方程(1.5.1)，设解的形式为 \(y_m = r^m\)（\(r\) 为常数），代入方程得：

\[r^{m+k} + \sum_{j=0}^{k-1} \alpha_j r^{m+j} = 0 \]

约去 \(r^m\)（\(r\neq0\)），得到特征方程：

\[\boldsymbol{\rho(r) = r^k + \sum_{j=0}^{k-1} \alpha_j r^j = 0} \tag{1.5.4} \]

其中 \(\rho(r)\) 称为第一特征多项式。

差分方程的通解

  • 若特征方程(1.5.4)有 \(k\) 个互不相同的根 \(r_1, r_2, \dots, r_k\)，则差分方程(1.5.1)的通解为：

    \[y_m = \sum_{i=1}^k C_i r_i^m \]

  • 若存在重根（如 \(r_1\) 是 \(s\) 重根），则通解中对应项为 \((C_1 + C_2 m + \dots + C_s m^{s-1}) r_1^m\)。

二、零稳定性（根条件）

2.1 零稳定性的定义

称线性多步法(1.4.1)是零稳定的，如果存在常数 \(C>0\) 和 \(h_0>0\)，使得对任意 \(0<h<h_0\)，方法的任意两个解 \(y_m\) 和 \(z_m\) 满足：

\[\max_{0\leq m\leq N} |y_m - z_m| \leq C \max_{0\leq j\leq k-1} |y_j - z_j| \]

其中 \(N = \frac{b-a}{h}\)。

核心含义：初始值的微小扰动，在整个计算过程中不会被无限放大，始终被初始扰动的常数倍控制。

2.2 零稳定性的充要条件（根条件）

定理1.5.1：线性多步法(1.4.1)是零稳定的，当且仅当它满足根条件：

特征方程 \(\rho(r)=0\) 的所有根的模都不超过1，且模为1的根都是单根。

根条件的直观解释

  • 若存在根 \(|r_i|>1\)：则 \(r_i^m\) 随 \(m\) 指数增长，初始扰动会被无限放大，方法零不稳定；
  • 若存在模为1的重根：则对应项为 \(m^s r_i^m\)，随 \(m\) 多项式增长，也会导致扰动无限放大，方法零不稳定；
  • 若所有根 \(|r_i|\leq1\)，且模为1的根为单根：则 \(r_i^m\) 有界，扰动不会无限放大，方法零稳定。

例子：Milne方法的零稳定性

Milne方法：\(y_{m+4} = y_m + \frac{4h}{3}(2f_{m+3}-f_{m+2}+2f_{m+1})\)

  • 特征方程：\(\rho(r) = r^4 - 1 = 0\)
  • 特征根：\(r=1, -1, i, -i\)，模均为1，且均为单根，满足根条件，零稳定。

例子：Simpson方法的零稳定性

Simpson方法：\(y_{m+2} = y_m + \frac{h}{3}(f_{m+2}+4f_{m+1}+f_m)\)

  • 特征方程：\(\rho(r) = r^2 - 1 = 0\)
  • 特征根：\(r=1, -1\)，模均为1，且均为单根，满足根条件，零稳定，但存在寄生根 \(r=-1\)，长区间计算会出现振荡。

三、局部截断误差与阶条件

3.1 局部截断误差的定义

定义1.5.1：设 \(y(x)\) 是初值问题的精确解，令 \(y_{m+j} = y(x_{m+j})\)（\(j=0,1,\dots,k-1\)），用方法(1.4.1)计算得到的 \(y_{m+k}\) 记为 \(y_{m+k}^*\)，则局部截断误差定义为：

\[\boldsymbol{R_{m,k} = y(x_{m+k}) - y_{m+k}^*} \]

若 \(R_{m,k} = O(h^{p+1})\)，则称方法是p阶精度的。

3.2 阶条件的等价形式

将精确解 \(y(x_{m+j})\) 代入方法(1.4.1)，并在 \(x_m\) 处做Taylor展开，整理得局部截断误差的展开式：

\[R_{m,k} = \sum_{n=0}^\infty \frac{h^n}{n!} C_n y^{(n)}(x_m) \]

其中系数 \(C_n\) 为：

\[\begin{cases} C_0 = 1 + \sum_{j=0}^{k-1} \alpha_j \\ C_n = \frac{1}{n!} \left( k^n + \sum_{j=0}^{k-1} \alpha_j j^n - n \sum_{j=0}^k \beta_j j^{n-1} \right), \quad n\geq1 \end{cases} \]

阶条件：要使方法为p阶精度，必须满足：

\[C_0 = C_1 = \dots = C_p = 0, \quad C_{p+1} \neq 0 \]

  • 若 \(C_0 = C_1 = 0\)，则方法与微分方程相容（\(p\geq1\)）。

3.3 线性多步法的阶数上限

定理1.5.10：

  • 显式k步线性多步法：最高可达到k阶精度；
  • 隐式k步线性多步法：最高可达到2k阶精度（如Gauss型隐式RK方法）。

四、收敛性：Lax等价定理

4.1 收敛性的定义

称线性多步法(1.4.1)是收敛的，如果当 \(h\to0\) 时，对任意初始值 \(y_j \to y(x_j)\)（\(j=0,1,\dots,k-1\)），都有：

\[\max_{0\leq m\leq N} |y(x_m) - y_m| \to 0 \]

4.2 收敛性的充要条件（Lax等价定理）

定理1.5.14：线性多步法(1.4.1)收敛的充要条件是：

  1. 相容：方法阶数 \(p\geq1\)（即 \(C_0 = C_1 = 0\)）；
  2. 零稳定：满足根条件。

核心推论

  • 若方法相容且零稳定，则整体截断误差为 \(O(h^p)\)，即p阶收敛；
  • 若方法不相容或零不稳定，则一定不收敛。

反例：Milne方法的收敛性

Milne方法：\(p=4\)（相容），且满足根条件（零稳定），因此收敛，但因存在寄生根 \(r=-1\)，长区间计算会出现数值振荡。

反例：一个不稳定的方法

考虑方法：\(y_{m+2} - 2y_{m+1} + 2y_m = h f_{m+1}\)

  • 特征方程：\(\rho(r) = r^2 - 2r + 2 = 0\)，根为 \(r=1\pm i\)，模为 \(\sqrt{2}>1\)，不满足根条件，零不稳定；
  • 尽管 \(C_0=0, C_1=0\)（相容），但方法不收敛，数值解会随步数指数增长。

五、绝对稳定性

5.1 绝对稳定性的定义

采用线性测试方程 \(y' = \lambda y\)（\(\text{Re}(\lambda)<0\)），将线性多步法代入得：

\[y_{m+k} + \sum_{j=0}^{k-1} \alpha_j y_{m+j} = z \sum_{j=0}^k \beta_j y_{m+j}, \quad z = \lambda h \]

定义稳定性多项式：

\[\boldsymbol{\pi(r,z) = \rho(r) - z \sigma(r) = 0} \]

其中 \(\sigma(r) = \sum_{j=0}^k \beta_j r^j\) 为第二特征多项式。

称方法在 \(z\) 点是绝对稳定的，如果 \(\pi(r,z)=0\) 的所有根的模都小于1；复平面上所有满足该条件的 \(z\) 的集合称为绝对稳定区域。

5.2 常用线性多步法的绝对稳定区间

方法 步数 阶数 实轴绝对稳定区间
显式Adams-Bashforth 1阶 1 1 \((-2, 0)\)
显式Adams-Bashforth 2阶 2 2 \((-1, 0)\)
显式Adams-Bashforth 4阶 4 4 \((-0.30, 0)\)
隐式Adams-Moulton 2阶 1 2 \((-\infty, 0)\)
隐式Adams-Moulton 4阶 3 4 \((-3, 0)\)
4阶Adams预测-校正 4 4 \((-0.7, 0)\)

5.3 绝对稳定性与零稳定性的区别

  • 零稳定性：\(h\to0\) 时的渐进稳定性，刻画初始扰动的长期行为；
  • 绝对稳定性：有限步长 \(h\) 下的稳定性，刻画误差在单步计算中的增长行为；
  • 零稳定是绝对稳定的前提，但零稳定的方法不一定绝对稳定（如Milne方法零稳定，但绝对稳定区间很小）。

六、核心知识点系统性归纳总结

表1 线性多步法核心理论结论对比

核心概念 核心定义 充要条件/关键结论
零稳定性 初始扰动不会随步数无限放大 充要条件：满足根条件（\(\rho(r)=0\) 的根模≤1，模1的根为单根）
相容性 方法与微分方程一致，局部截断误差 \(O(h^2)\) 充要条件：\(C_0 = C_1 = 0\)（\(p\geq1\)）
收敛性 \(h\to0\) 时数值解收敛于精确解 充要条件：相容 + 零稳定（Lax等价定理）
局部截断误差 前k步精确时，单步计算的误差 \(R_{m,k} = O(h^{p+1})\) ⇒ 方法为p阶精度
绝对稳定性 有限步长下误差不会指数增长 充要条件：\(\pi(r,z)=0\) 的所有根模<1
阶数上限 方法可达到的最高精度阶数 显式k步法：最高k阶；隐式k步法：最高2k阶

表2 常用线性多步法的稳定性与收敛性

方法名称 步数 阶数 特征方程 根条件满足？ 零稳定？ 收敛？ 实轴绝对稳定区间
显式Euler 1 1 \(r-1=0\) 是 是 是 \((-2,0)\)
2步2阶Adams-Bashforth 2 2 \(r^2 - r = 0\) 是 是 是 \((-1,0)\)
4步4阶Adams-Bashforth 4 4 \(r^4 - r^3 = 0\) 是 是 是 \((-0.30,0)\)
梯形法 1 2 \(r-1=0\) 是 是 是 \((-\infty,0)\)
3步4阶Adams-Moulton 3 4 \(r^3 - r^2 = 0\) 是 是 是 \((-3,0)\)
Milne方法 4 4 \(r^4 - 1 = 0\) 是 是 是 不稳定区间大
Simpson方法 2 4 \(r^2 - 1 = 0\) 是 是 是 不稳定区间大
不稳定方法 2 1 \(r^2 - 2r + 2 = 0\) 否 否 否 无

表3 线性多步法与单步法的核心理论对比

特性 线性多步法 单步法（Runge-Kutta）
自启动 不可自启动，需单步法提供初始值 可自启动，仅需初始值 \(y_0\)
零稳定性 需满足根条件，存在零不稳定的方法 所有单步法均零稳定
收敛性 相容+零稳定 ⇔ 收敛 相容 ⇔ 收敛（单步法天然零稳定）
绝对稳定性 隐式方法绝对稳定区间大，显式方法小 显式RK绝对稳定区间有限，隐式RK可A-稳定
阶数上限 显式k步法k阶，隐式k步法2k阶 显式s级RK最多s阶
计算效率 同阶精度下，每步仅算1次f，效率极高 同阶精度下，每步需多次算f，效率较低

七、关键定理与结论速记

  1. Lax等价定理：线性多步法收敛 ⇔ 相容 + 零稳定。
  2. 根条件：零稳定 ⇔ \(\rho(r)=0\) 的根模≤1，模1的根为单根。
  3. 阶条件：p阶精度 ⇔ \(C_0 = \dots = C_p = 0\)，\(C_{p+1} \neq 0\)。
  4. 绝对稳定：\(\pi(r,z)=0\) 的所有根模<1 ⇔ 方法在z点绝对稳定。
  5. 阶数上限：显式k步法≤k阶，隐式k步法≤2k阶。

posted on 2026-03-22 08:17  Indian_Mysore  阅读(77)  评论(0)    收藏  举报

导航