常微分方程初值问题
许多物理、工程和经济学中的动态系统都可以用常微分方程 (Ordinary Differential Equations, ODEs) 来描述。一个典型的一阶ODE具有以下形式：
$$y'(t) = f(t, y(t)) $$
这个方程描述了变量 $y$ 随时间 $t$ 变化的速率。
然而，仅有这个方程不足以确定系统的具体行为，因为它有无穷多个可能的解。为了得到一个唯一的解，我们还需要一个初始状态，即初值条件 (Initial Condition)：
$$y(t_0) = y_0 $$
将ODE和初值条件结合在一起，就构成了一个初值问题 (Initial Value Problem, IVP)。我们的目标就是找到函数 $y(t)$，它在 $t \ge t_0$ 的区间上满足上述两个条件。
解的存在性与唯一性理论
在尝试用数值方法求解之前，我们必须先确定一个问题是否有唯一的、行为良好的解。
Lipschitz 连续性
Lipschitz连续性是保证ODE解行为良好的一个关键概念。
定义：如果存在一个正常数 $L$（称为Lipschitz常数），使得对于任意的 $y_1, y_2$，函数 $f(t, y)$ 都满足以下不等式，则称函数 $f(t, y)$ 关于其第二个变量 $y$ 是Lipschitz连续的：
$$|f(t, y_1) - f(t, y_2)| \le L |y_1 - y_2| $$
直观上，这个条件限制了函数 $f$ 随 $y$ 变化的“陡峭”程度。它保证了当 $y$ 发生微小变化时，$y$ 的变化率 $f(t, y)$ 不会发生爆炸性增长。
Cauchy-Lipschitz 定理 (存在唯一性定理)
这个定理为IVP解的存在性和唯一性提供了坚实的理论基础。
定理：对于初值问题 $y'(t) = f(t, y(t)), y(t_0) = y_0$，如果函数 $f(t, y)$ 在一个包含点 $(t_0, y_0)$ 的区域内是关于 $t$ 连续的，且关于 $y$ 是Lipschitz连续的，那么在 $t_0$ 附近的一个时间区间内，该IVP存在一个唯一的、连续可微的解 $y(t)$。
这个定理告诉我们，只要函数 $f$ 足够“温和”，我们所求解的目标就是明确且唯一的。
离散化：从连续到离散的桥梁
数值方法的核心思想是用一系列离散点上的近似值来代替连续的真实解。我们设定一系列时间点 $t_0, t_1, t_2, \ldots$，其中步长为 $h = t_{n+1} - t_n$，并计算在这些点上的近似解 $y_0, y_1, y_2, \ldots$，使得 $y_n \approx y(t_n)$。
这个过程的理论基础来源于牛顿-莱布尼茨公式。将ODE $y'(t) = f(t, y(t))$ 在区间 $[t_n, t_{n+1}]$ 上积分，我们得到：
$$\int_{t_n}^{t_{n+1}} y'(t) dt = y(t_{n+1}) - y(t_n) $$
因此，真实解满足：
$$y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t, y(t)) dt $$
这个公式是精确的。数值方法的本质区别，就在于如何近似右侧的积分项 $\int_{t_n}^{t_{n+1}} f(t, y(t)) dt$。
经典的数值求解方法
前向欧拉法 (Forward Euler)
最简单的近似方法是用积分区间左端点的值来近似整个积分，即 $\int_{t_n}^{t_{n+1}} f(t, y(t)) dt \approx h \cdot f(t_n, y(t_n))$。
由此得到前向欧拉法的迭代公式：
$$y_{n+1} = y_n + h f(t_n, y_n) $$
这是一个显式方法，因为 $y_{n+1}$ 的计算直接由已知的 $y_n$ 给出。
后向欧拉法 (Backward Euler)
如果我们用积分区间右端点的值来近似积分，即 $\int_{t_n}^{t_{n+1}} f(t, y(t)) dt \approx h \cdot f(t_{n+1}, y(t_{n+1}))$。
迭代公式变为：
$$y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) $$
这是一个隐式方法，因为待求的 $y_{n+1}$ 同时出现在方程的两边。求解它通常需要使用牛顿法等非线性方程求解器。
Crank-Nicolson 方法 (梯形法)
为了获得更高的精度，我们可以使用梯形法则来近似积分，即取左右端点的平均值：
$$y_{n+1} = y_n + \frac{h}{2} [f(t_n, y_n) + f(t_{n+1}, y_{n+1})] $$
Crank-Nicolson方法同样是隐式的，但其精度通常高于前向和后向欧拉法。
欧拉-柯西 预估-校正法 (Predictor-Corrector)
这个方法试图在不进行复杂隐式求解的情况下提高精度。它分为两步：
- 预估 (Predict)：使用前向欧拉法得到一个对 $y_{n+1}$ 的初步估计值 $\tilde{y}_{n+1}$。
$$\tilde{y}_{n+1} = y_n + h f(t_n, y_n) $$
- 校正 (Correct)：使用这个估计值来近似 $f$ 在 $t_{n+1}$ 处的值，并用梯形法则进行校正。
$$y_{n+1} = y_n + \frac{h}{2} [f(t_n, y_n) + f(t_{n+1}, \tilde{y}_{n+1})] $$
这种方法是显式的，但其精度与隐式的Crank-Nicolson法相当。
单步多阶段方法：Runge–Kutta（RK）
Runge–Kutta 属于单步方法（只依赖最近一步的 $(t_n,y_n)$），但在同一步内通过多个“阶段”评估 $f(t,y)$ 来提高阶数与稳定性。
经典二阶RK（举例）
-
Heun（显式梯形）法：
- $k_1 = f(t_n, y_n)$
- $k_2 = f(t_n+h,\, y_n + h k_1)$
- $y_{n+1} = y_n + \tfrac{h}{2}(k_1 + k_2)$
-
中点法：
- $k_1 = f(t_n, y_n)$
- $k_2 = f(t_n+\tfrac{h}{2},\, y_n + \tfrac{h}{2} k_1)$
- $y_{n+1} = y_n + h\,k_2$
两者均为二阶方法：局部截断误差 $O(h^{3})$，全局误差 $O(h^{2})$。
经典四阶RK（RK4）
$$\begin{aligned} k_1 &= f(t_n, y_n),\\ k_2 &= f\bigl(t_n+\tfrac{h}{2},\, y_n + \tfrac{h}{2} k_1\bigr),\\ k_3 &= f\bigl(t_n+\tfrac{h}{2},\, y_n + \tfrac{h}{2} k_2\bigr),\\ k_4 &= f\bigl(t_n+h,\, y_n + h k_3\bigr),\\ y_{n+1} &= y_n + \tfrac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{aligned} $$
RK4 的局部截断误差为 $O(h^{5})$、全局误差为 $O(h^{4})$，在非刚性问题中是极其常用的高性价比方法。
Butcher 表与阶条件（简述）
RK 方法可用 Butcher 表表征各阶段系数。阶条件由代数等式约束这些系数，使方法达到既定阶数（此处略）。
自适应步长（嵌入式对）
常用的 RKF45、Dormand–Prince（DOPRI5）等提供两个相近阶数的解 $y^{[p]}$ 与 $y^{[p-1]}$ 以估计单步误差 $\Delta$：
$$\Delta = \lVert y^{[p]} - y^{[p-1]} \rVert. $$
据此调整步长，例如
$$h_{\text{new}} = h\;\mathrm{safety}\;\left(\frac{\mathrm{tol}}{\Delta}\right)^{\frac{1}{p+1}}, $$
并约束放大/缩小因子在合理区间内，以保证效率与稳健性。
线性测试方程与稳定性区域
考虑 $y' = \lambda y$，令 $z = h\lambda$，RK 方法对应的放大因子 $R(z)$ 决定了稳定性区域 $\{z:\,|R(z)|\le 1\}$。显式 RK 的稳定域有界，步长需满足 $h$ 足够小；刚性问题宜选隐式 RK（如 Gauss–Legendre）或下述 BDF。
线性多步方法（LMM）
线性多步法在一步中使用多个历史点来提高阶数，一般形式：
$$\sum_{j=0}^{k} \alpha_j\, y_{n-j} = h \sum_{j=0}^{k} \beta_j\, f(t_{n-j}, y_{n-j}). $$
一致性与零稳定性（根条件）
- 一致性（必要条件示例）：$\sum_{j=0}^k \alpha_j = 0$，且 $\sum_{j=0}^k j\alpha_j = \sum_{j=0}^k \beta_j$。
- 零稳定性：特征多项式 $\rho(\xi)=\sum_{j=0}^k \alpha_j\xi^{k-j}$ 的根均在单位圆内，且在单位圆上的根为单根。零稳定性 + 一致性 ⇒ 收敛。
Adams 家族（积分型）
-
显式 Adams–Bashforth（AB）二阶：
$$y_{n+1} = y_n + h\left(\tfrac{3}{2} f_n - \tfrac{1}{2} f_{n-1}\right). $$
-
隐式 Adams–Moulton（AM）二阶（梯形/Crank–Nicolson）：
$$y_{n+1} = y_n + \tfrac{h}{2}\bigl(f_{n+1} + f_n\bigr). $$
AB/AM 可组成预测–校正对：先用 AB 预测 $\tilde y_{n+1}$，再用 AM 校正 $y_{n+1}$，既保持较低成本又提升精度与稳定性。
后向差分公式（BDF，微分型）
-
BDF1（即后向欧拉，L-稳定）：
$$\frac{y_{n+1} - y_n}{h} = f(t_{n+1}, y_{n+1}). $$
-
BDF2（A-稳定）：
$$\frac{3y_{n+1} - 4y_n + y_{n-1}}{2h} = f(t_{n+1}, y_{n+1}). $$
BDF 对刚性问题表现优异。一般地，BDF 阶数超过 2 时不再 A-稳定，但仍常用于刚性问题的实际工程求解。
启动策略与实现要点
- 多步法需要前 $k$ 个起始值，通常用高阶 RK 在初段“引导”。
- 隐式 AM/BDF 需解非线性方程，可用牛顿或简化牛顿（保持雅可比）。
- 自适应步长/阶数：Adams 与 BDF 系列普遍支持按误差估计调整 $h$ 与阶数。
选择建议（非刚性 vs 刚性）
- 平滑、非刚性：显式 RK（如 RK4、DOPRI5）配合自适应步长，简单稳健。
- 刚性：隐式阶梯（CN）、BDF、隐式 RK；必要时使用稀疏/预条件线性代数加速。
收敛性、稳定性与误差
一个好的数值方法必须是收敛的，即当步长 $h \to 0$ 时，数值解要趋近于真实解。收敛性由两个关键因素决定：一致性和稳定性。
误差来源
- 局部截断误差 (Local Truncation Error)：假设 $y_n$ 是精确的，即 $y_n = y(t_n)$，在计算一步之后，数值解 $y_{n+1}$ 与真实解 $y(t_{n+1})$ 之间的误差。它衡量了方法本身的近似程度。
- 全局误差 (Global Error)：在经过多步计算后，由于每一步的局部误差不断累积，最终在 $t_n$ 时刻的数值解 $y_n$ 与真实解 $y(t_n)$ 之间的总误差。
一致性 (Consistency)
一致性衡量的是局部截断误差的大小。如果一个方法的局部截断误差为 $O(h^{p+1})$，我们称这个方法是 $p$ 阶一致的。
- 定义：当步长 $h \to 0$ 时，如果局部截断误差也趋于0，则称该方法是一致的。
- 可靠性：一致性保证了在单步计算中，我们的方法是对微分方程的“忠实”近似。
稳定性 (Stability)
稳定性（或称零稳定性）衡量的是误差的累积行为。
- 定义：如果一个方法在计算过程中不会将微小的扰动（如舍入误差或局部误差）无限放大，导致最终结果发散，那么该方法是稳定的。
- 重要性：不稳定的方法即使在理论上是一致的，在实际计算中也毫无用处，因为误差会淹没真实解。
收敛定理
数值分析中的一个基本定理指出：
$$\Large \boxed{ \textbf{一致性}\; \boldsymbol{+}\; \textbf{稳定性}\; \boldsymbol{\Rightarrow}\; \textbf{收敛} } $$
这意味着，一个数值方法要想收敛，它必须既能忠实地近似原方程（一致性），又能有效地控制误差的传播（稳定性）。对于我们上面介绍的几种方法，它们都是收敛的。