📑 目录

1. 微分方程数值解概览
2. Euler 法与 Taylor 展开法
3. Runge-Kutta 方法族
4. Predictor-Corrector 多步法
5. 微分方程组与 Lorenz 混沌
6. 高阶微分方程的降阶处理
7. 刚性微分方程（Stiff ODE）
8. 边值问题与打靶法（Shooting）
9. Numerov 方法
10. Schrödinger 方程数值求解
11. Green 函数与 Wronskian
12. 有限差分法（Finite Difference）
13. 应用实例：TOV 方程与中子星
14. 总结对比表

1. 微分方程数值解概览

两大问题类型

💡 物理直觉：初值问题是"知道现在，推演未来"（动力学），边值问题是"知道两端，反推中间"（稳态分布）。

2. Euler 法与 Taylor 展开法

2.1 Euler 法 — 最朴素的一阶方法

核心思想：用切线近似曲线。

对于一阶 ODE：

等间距离散化：取步长 $h$，定义 $t_i = t_0 + i \cdot h$，$i = 0,1,2,...,N$

Euler 递推公式：

说人话：下一步 = 当前值 + 步长 × 当前斜率。就像盲人走路，每步朝着当前方向走固定距离，完全不看前面弯不弯。

📐 数值例子： $y' = -2y + 1,\ y(0)=1$

精确解： $y(t) = 0.5 + 0.5e^{-2t}$

# Euler 法手算演示，h=0.2
import numpy as np

def f(t, y): return -2*y + 1
y0, h = 1.0, 0.2
t_vals = np.arange(0, 1.01, h)
y_euler = [y0]

for i in range(len(t_vals)-1):
    y_next = y_euler[-1] + h * f(t_vals[i], y_euler[-1])
    y_euler.append(y_next)

观察：Euler 法在 $t=1.0$ 时误差约 12%。步长太大！

2.2 高阶 Taylor 展开法

递推公式：

局部截断误差 (Local Truncation Error)：

• Euler 法（1阶 Taylor）：$O(h^2)$ 每步
• $k$ 阶 Taylor：$O(h^{k+1})$ 每步

对 $y'=-2y+1$ 用二阶 Taylor

$y'' = -2y' = -2(-2y+1) = 4y-2$

# 二阶 Taylor，同一步长 h=0.2
y_taylor2 = [y0]
for i in range(len(t_vals)-1):
    y = y_taylor2[-1]
    yp = -2*y + 1          # y'
    ypp = 4*y - 2          # y''
    y_next = y + h*yp + h**2/2 * ypp
    y_taylor2.append(y_next)

# t=1.0: y≈0.559592 → 误差从 12% 降到约 1.4%！

💡 结论：提高阶数比缩小步长更有效！二阶 Taylor 用 $h=0.2$ 的精度已经超过 Euler 用 $h=0.1$。

3. Runge-Kutta 方法族

在科学计算、航空航天轨道计算、乃至游戏物理引擎中，90% 以上的常微分方程初值问题，其底层默认的算法都是龙格-库塔法。

3.1 核心思想

Taylor 展开法需要手动求各阶导数

（麻烦！ + 阶数越高 手动求偏导的项数呈现指数级爆炸 —— 工程上不可以承受的 + 导数无具体解析解 需要数值手段）

Runge-Kutta 的巧思：

能否不求导，仅仅是通过在当前步长的不同位置 多算几次斜率，然后将其“加权平均”，来拼凑出来和Taylor一样的高阶精度效果

用多点函数值的组合来模拟高阶导数效果，无需显式求导。 我们无需去求那些复杂的偏导数，我们只需要在区间 $[t_{i}, t_{i+1}]$ 之间 ，多选取

一般形式：

其中 k就是我们在这一步里，通过不同位置的“试探”得到的斜率。每一个k，都依赖于前面已经算出来的斜率

这些系数 $a_{jl}$、 $b_j$ 和 $c_j$ 是通过与泰勒级数展开进行复杂的代数比对（Taylor matching）强行解出来

3.2 二阶 Runge-Kutta：中点法与改进 Euler 法

🔹 中点法 (Midpoint Method) 先探路 再迈步

说人话：先用半步走到的中点斜率，再拿这个斜率走一整步。比直接走精度高了一阶。

k1就是 起点的斜率 ， k2就是中点处预测的斜率

你蒙着眼睛在一座山上往下走。如果直接用起点斜率 $k_1$ 走一整步（欧拉法），由于山坡在弯曲，你很容易偏离轨道。 中点法的策略是：我先用起点的斜率 $k_1$ 试探着走半步（ $\frac{h}{2}$），到达“中点”。我测出中点处的斜率 $k_2$，然后我回到起点，手持这个中点斜率 $k_2$ 跨出完整的一整步。由于中点斜率代表了这一小段路程的“平均趋势”，这样跨出的一步，精度直接从一阶跃升到了二阶（全局误差 $O(h^2)$）！

🔹 改进 Euler 法 (Modified Euler / Heun's Method)

说人话：取"起点斜率"和"用 Euler 走到终点的斜率"的平均。相当于预测-校正。

中点法 是去中间探路 而 改进的欧拉法 则是两头兼顾

可以先用 起点的斜率 k1 算一个“预期的终点”，接着 我们测量这个预期终点处斜率k2，最后将起点与终点的斜率进行算术平均 得出最终迈步的依据 —— 预测-校正

Predictor- Corrector 思想！！！

3.3 四阶 Runge-Kutta (RK4) — 工业标准

公式（经典 RK4）：

$\begin{aligned} k_1 &= f(t_i, y_i) \\ \space \\ k_2 &= f\left(t_i + \frac{h}{2},\ y_i + \frac{h}{2}k_1\right) \\ \space \\ k_3 &= f\left(t_i + \frac{h}{2},\ y_i + \frac{h}{2}k_2\right) \\ \space \\ k_4 &= f(t_i + h,\ y_i + h k_3) \\ \space \\ y_{i+1} &= y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}$

本质上其代数逻辑 与 Simpson 积分 惊人的相似！！！

与数值积分里的 Simpson 公式权重（ $\frac{1}{6}, \frac{4}{6}, \frac{1}{6}$）极其相似！

$k_1$（起点斜率）权重为 $1/6$；

$k_2, k_3$（两个在中点不同试探方式得到的斜率）权重合起来是 $4/6$；

$k_4$（终点斜率）权重为 $1/6$。

本质

是常微分方程数值解和数值积分 在深层代数逻辑上 合流的必然结果

一、 终极源头：将 ODE 降维成“纯积分”

要理解为什么它们会发生碰撞，我们先做一个数学实验： 假设我们有一个一阶常微分方程： $\frac{dy}{dt} = f(t, y), \quad y(t_i) = y_i$

如果我们对两边在区间 $[t_i, t_i + h]$ 上进行精确求导积分，由微积分基本定理可以得到： $y(t_i + h) - y(t_i) = \int_{t_i}^{t_i+h} f(t, y(t)) dt$

所以，下一步的精确值一定是： $y_{i+1} = y_i + \int_{t_i}^{t_i+h} f(t, y(t)) dt$

这是一个无法直接计算的积分，因为被积函数里的 $y(t)$ 本身就是我们正在求解的未知数！

降维假设：如果斜率不依赖于 $y$

现在，我们假设一种最简单、最特殊的物理场景：斜率 $f(t, y)$ 仅仅与时间 $t$ 有关，与 $y$ 无关。 比如，一个物体的加速度只随时间变化： $\frac{dy}{dt} = f(t)$。

此时，上面的方程直接变成了一个纯积分问题： $y_{i+1} = y_i + \int_{t_i}^{t_i+h} f(t) dt$

二、 辛普森公式（Simpson's Rule）是如何处理这个积分的？

如果要计算 $\int_{t_i}^{t_i+h} f(t) dt$，辛普森公式的思路是：在区间 $[t_i, t_i+h]$ 里，取三个等距点来做二次抛物线插值：

1. 左端点： $t_i$，对应的函数值为 $f(t_i)$；
2. 中点： $t_i + \frac{h}{2}$，对应的函数值为 $f(t_i + \frac{h}{2})$；
3. 右端点： $t_i + h$，对应的函数值为 $f(t_i + h)$。

根据辛普森公式的求积权重（$1:4:1$ 比例），这个积分的近似值为： $\int_{t_i}^{t_i+h} f(t) dt \approx \frac{h}{6} \left[ f(t_i) + 4 f\left(t_i + \frac{h}{2}\right) + f(t_i + h) \right]$

因此，我们得到下一步的预测值为：

$y_{i+1} \approx y_i + \frac{h}{6} \left[ f(t_i) + 4 f\left(t_i + \frac{h}{2}\right) + f(t_i + h) \right]$

三、 经典 RK4 此时在干什么？

现在，我们把相同的“降维假设”（$f$ 不依赖于 $y$）代入到经典 RK4 的公式里，看看它会退化成什么样3：

1. 计算 $k_1$（起点斜率）： $k_1 = f(t_i, y_i) \implies k_1 = f(t_i)$
2. 计算 $k_2$（中点试探斜率 1）： $k_2 = f\left(t_i + \frac{h}{2},\ y_i + \frac{h}{2}k_1\right) \implies k_2 = f\left(t_i + \frac{h}{2}\right)$
3. 计算 $k_3$（中点试探斜率 2）： $k_3 = f\left(t_i + \frac{h}{2},\ y_i + \frac{h}{2}k_2\right) \implies k_3 = f\left(t_i + \frac{h}{2}\right)$
4. 计算 $k_4$（终点试探斜率）： $k_4 = f(t_i + h,\ y_i + h k_3) \implies k_4 = f(t_i + h)$

因为 $f$ 与 $y$ 无关，我们惊奇地发现：在中点算出来的两个试探斜率 $k_2$ 和 $k_3$ 居然完全相等了2！ $k_2 = k_3 = f\left(t_i + \frac{h}{2}\right)$

现在，我们把它们代入 RK4 的最后一步加权平均式中： $\begin{aligned} y_{i+1} &= y_i + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \\ &= y_i + \frac{h}{6} \left[ f(t_i) + 2f\left(t_i + \frac{h}{2}\right) + 2f\left(t_i + \frac{h}{2}\right) + f(t_i + h) \right] \\ &= y_i + \frac{h}{6} \left[ f(t_i) + 4f\left(t_i + \frac{h}{2}\right) + f(t_i + h) \right] \end{aligned}$

两套公式完全合二为一，分毫不差！

四、 这个“巧合”背后的深刻含义

1. RK4 是 辛普森公式在 ODE 领域的“非线性推广”

• 辛普森公式：只能处理函数只与时间 $t$ 相关的静态积分。

• RK4：当函数不仅与 $t$ 有关，还与当前状态 $y$ 动态耦合时，RK4 通过巧妙的多步迭代预测（用 $k_1$ 预测出 $y_{i} + \frac{h}{2}k_1$，进而算 $k_2$；再用 $k_2$ 修正出 $y_{i} + \frac{h}{2}k_2$，进而算 $k_3$），在非线性空间里重建了辛普森公式所需要的“中点斜率”（在泰勒展开的代数精度上，RK4 巧妙地通过这两步试探，将非线性耦合带来的高阶误差一项一项全部抵消掉，从而在代数精度上达到了与辛普森公式完全匹配的四阶精度 $O(h^4)$）

• 因此，RK4 本质上就是“带有动态状态修正的辛普森积分”。

补充

一、 核心痛点：辛普森公式需要“真正的中点斜率”

如果我们把 ODE 的精确积分形式写出来： $y_{i+1} = y_i + \int_{t_i}^{t_i+h} f(t, y(t)) dt$ 如果我们想用辛普森公式来算这个积分，我们需要知道三个值：

1. 起点斜率： $f(t_i, y_i)$ —— 已知，就是 $k_1$。
2. 终点斜率： $f(t_i + h, y(t_i+h))$ —— 未知，因为我们不知道终点的精确状态 $y(t_i+h)$。
3. 中点真正的斜率： $f\left(t_i + \frac{h}{2}, y\left(t_i + \frac{h}{2}\right)\right)$ —— 未知，因为我们不知道中点精确状态 $y\left(t_i + \frac{h}{2}\right)$。

因为我们不知道中点精确状态 $y\left(t_i + \frac{h}{2}\right)$，所以我们无法直接算出真正的中点斜率。

二、 RK4 是如何用 $k_1, k_2, k_3, k_4$ 逐步逼近这个“真中点斜率”的？

RK4 的精妙之处在于，它通过三次接力（预测-校正-再预测），把这个“真正的中点斜率”给强行“套”了出来。

根据泰勒展开，函数 $y(t)$ 在 $t_i$ 处的半步（$\frac{h}{2}$）展开为： $y\left(t_i + \frac{h}{2}\right) = y(t_i) + \left(\frac{h}{2}\right) y'(t_i) + \frac{1}{2!} \left(\frac{h}{2}\right)^2 y''(t_i) + O(h^3)$ 整理一下系数： $y\left(t_i + \frac{h}{2}\right) = y_i + \frac{h}{2} y'(t_i) + \frac{h^2}{8} y''(t_i) + O(h^3) \quad \text{--- (式 A)}$

因为我们的常微分方程是 $y' = f(t, y)$，所以利用全微分（链式法则）求二阶导数 $y''$： $y''(t) = \frac{d}{dt}[f(t, y(t))] = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial y} \cdot \frac{dy}{dt} = f_t + f_y f$ （为了书写简洁，我们用偏导数简写 $f_t = \frac{\partial f}{\partial t}$，$fy = \frac{\partial f}{\partial y}$，且它们都在 $(ti, y_i)$ 处取值）。

把 $y'(t_i) = f$ 和 $y''(t_i) = f_t + f_y f$ 代回 (式 A)，得到真正的中点状态的精确展开式： $y\left(t_i + \frac{h}{2}\right) = y_i + \frac{h}{2} f + \frac{h^2}{8} (f_t + f_y f) + O(h^3) \quad \text{--- (式 A-精确版)}$

在 RK4 中，我们第一步用 $k_1$（即 $f$）预测的中点状态为： $y_{\text{mid}, 1} = y_i + \frac{h}{2} k_1 = y_i + \frac{h}{2} f$

我们用它去计算 $k_2$： $k_2 = f\left(t_i + \frac{h}{2}, \ y_i + \frac{h}{2} f\right)$

这是一个二元函数 $f(t, y)$。我们将它在点 $(t_i, y_i)$ 处进行二元泰勒展开。 二元泰勒展开公式为： $f(t_i + \Delta t, \ y_i + \Delta y) = f(t_i, y_i) + \Delta t \cdot f_t + \Delta y \cdot f_y + O(\Delta^2)$

在 $k_2$ 的计算中：

• 时间增量 $\Delta t = \frac{h}{2}$
• 状态增量 $\Delta y = \frac{h}{2} f$

代入二元泰勒展开式： $k_2 = f + \left(\frac{h}{2}\right) f_t + \left(\frac{h}{2} f\right) f_y + O(h^2)$ 提取公因子 $\frac{h}{2}$ $k_2 = f + \frac{h}{2} (f_t + f_y f) + O(h^2) \quad \text{--- (式 B)}$

现在，我们用这个含有更高阶误差信息的 $k_2$ 重新代入到第二次中点预测值 $y_{\text{mid}, 2}$ 中： $y_{\text{mid}, 2} = y_i + \frac{h}{2} k_2$

把 (式 B) 的展开式直接代入到 $k_2$ 的位置： $y_{\text{mid}, 2} = y_i + \frac{h}{2} \left[ f + \frac{h}{2} (f_t + f_y f) + O(h^2) \right]$

把括号外面的 $\frac{h}{2}$ 乘进去： $y_{\text{mid}, 2} = y_i + \frac{h}{2} f + \frac{h^2}{4} (f_t + f_y f) + O(h^3) \quad \text{--- (式 C)}$

在 RK4 公式中，第三步是在第二次中点预测状态 $y_{\text{mid}, 2}$ 处计算斜率： $k_3 = f\left(t_i + \frac{h}{2}, \ y_{\text{mid}, 2}\right)$

我们将 (式 C) 的结果： $y_{\text{mid}, 2} = y_i + \frac{h}{2} f + \frac{h^2}{4} (f_t + f_y f) + O(h^3)$ 代入到 $k_3$ 中： $k_3 = f\left(t_i + \frac{h}{2}, \ y_i + \frac{h}{2} f + \frac{h^2}{4} (f_t + f_y f) + O(h^3)\right)$

此时，我们需要对二元函数 $f(t_i + \Delta t, \ y_i + \Delta y)$ 进行泰勒展开，其中：

• 时间增量 $\Delta t = \frac{h}{2}$
• 状态增量 $\Delta y = \frac{h}{2} f + \frac{h^2}{4} (f_t + f_y f) + O(h^3)$

为了计算到 $O(h^2)$ 项，我们的二元泰勒展开需要保留到二阶偏导数项： $f(t_i + \Delta t, \ y_i + \Delta y) = f + \Delta t f_t + \Delta y f_y + \frac{1}{2} \left[ \Delta t^2 f_{tt} + 2 \Delta t \Delta y f_{ty} + \Delta y^2 f_{yy} \right] + O(h^3)$ （注：这里所有的偏导数 $f_t, f_y, f_{tt}, f_{ty}, f_{yy}$ 都在 $(t_i, y_i)$ 处取值）。

我们分项把它们代入展开：

• 一阶时间项： $\Delta t f_t = \frac{h}{2} f_t$
• 一阶状态项： $\Delta y f_y = \left[ \frac{h}{2} f + \frac{h^2}{4} (f_t + f_y f) + O(h^3) \right] f_y = \frac{h}{2} f f_y + \frac{h^2}{4} f_y (f_t + f_y f) + O(h^3)$

• 二阶偏导数项（由于我们只需要保留到 $O(h^2)$，所以 $\Delta y$ 里的 $O(h^2)$ 项在平方后会变成 $O(h^3)$ 以上，可以省去）：

  $\frac{1}{2} \Delta t^2 f_{tt} = \frac{1}{2} \left(\frac{h}{2}\right)^2 f_{tt} = \frac{h^2}{8} f_{tt}$ $\frac{1}{2} \cdot 2 \Delta t \Delta y f_{ty} = \left(\frac{h}{2}\right) \left(\frac{h}{2} f + O(h^2)\right) f_{ty} = \frac{h^2}{4} f f_{ty}$ $\frac{1}{2} \Delta y^2 f_{yy} = \frac{1}{2} \left(\frac{h}{2} f + O(h^2)\right)^2 f_{yy} = \frac{h^2}{8} f^2 f_{yy}$

在 RK4 中，我们用最新的精细中点斜率 $k_3$ 跨出完整的一整步，预测终点的状态 $y_{\text{end}}$： $y_{\text{end}} = y_i + h k_3$

将 (式 D) 代入 $y_{\text{end}}$ 的计算中： $y_{\text{end}} = y_i + h \left( f + \frac{h}{2} (f_t + f_y f) + O(h^2) \right)$ 把 $h$ 乘进去，得到预测的终点状态： $y_{\text{end}} = y_i + h f + \frac{h^2}{2} (f_t + f_y f) + O(h^3) \quad \text{--- (式 E)}$

现在，我们在终点计算最后一个斜率 $k_4$： $k_4 = f(t_i + h, \ y_{\text{end}}) = f\left(t_i + h, \ y_i + h f + \frac{h^2}{2} (f_t + f_y f) + O(h^3)\right)$

我们再次对二元函数进行展开，此时：

• 时间增量 $\Delta t = h$
• 状态增量 $\Delta y = h f + \frac{h^2}{2} (f_t + f_y f) + O(h^3)$

代入二元泰勒展开公式： $f(t_i + \Delta t, \ y_i + \Delta y) = f + \Delta t f_t + \Delta y f_y + \frac{1}{2} \left[ \Delta t^2 f_{tt} + 2 \Delta t \Delta y f_{ty} + \Delta y^2 f_{yy} \right] + O(h^3)$

• 一阶时间项： $\Delta t f_t = h f_t$
• 一阶状态项： $\Delta y f_y = \left[ h f + \frac{h^2}{2} (f_t + f_y f) + O(h^3) \right] f_y = h f f_y + \frac{h^2}{2} f_y (f_t + f_y f) + O(h^3)$
• 二阶偏导数项（省去 $O(h^3)$ 阶）： $\frac{1}{2} \Delta t^2 f_{tt} = \frac{1}{2} h^2 f_{tt}$ $\frac{1}{2} \cdot 2 \Delta t \Delta y f_{ty} = h (h f + O(h^2)) f_{ty} = h^2 f f_{ty}$ $\frac{1}{2} \Delta y^2 f_{yy} = \frac{1}{2} (h f + O(h^2))^2 f_{yy} = \frac{1}{2} h^2 f^2 f_{yy}$

• $h$ 的一阶项： $h f_t + h f f_y = h (f_t + f_y f)$
• $h$ 的二阶项： $\frac{h^2}{2} f_y (f_t + f_y f) + \frac{h^2}{2} f_{tt} + h^2 f f_{ty} + \frac{h^2}{2} f^2 f_{yy}$ 我们通分并提取公因子 $\frac{h^2}{2}$： $\frac{h^2}{2} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right]$

因此，我们得到了 $k_4$ 的完整展开式： $k_4 = f + h (f_t + f_y f) + \frac{h^2}{2} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right] + O(h^3) \quad \text{--- (式 F)}$

现在，四个斜率 $k_1$、$k2$、$k3$、$k_4$ 已经全部列阵完毕。我们把它们毫无保留地写在一起：

1. $k_1 = f$
2. $k_2 = f + \frac{h}{2} (f_t + f_y f) + \frac{h^2}{8} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} \right] + O(h^3)$
3. $k_3 = f + \frac{h}{2} (f_t + f_y f) + \frac{h^2}{8} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + 2 f_y (f_t + f_y f) \right] + O(h^3)$
4. $k_4 = f + h (f_t + f_y f) + \frac{h^2}{2} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right] + O(h^3)$

现在，我们将它们代入 RK4 最终的加权表达式： $y_{i+1} = y_i + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4)$

为了让你看清楚每一部分是如何抵消和合并的，我们把求和部分 $k_1 + 2k_2 + 2k_3 + k_4$ 按 $h$ 的阶数 拆成三部分算：

1. 零阶项（常数 $f$ 的加权）

$\text{零阶项和} = f + 2(f) + 2(f) + f = 6f$ 代回最终公式中，再乘以外面的 $\frac{h}{6}$： $\frac{h}{6} (6f) = \mathbf{h f}$

2. 一阶项（含 $f_t + f_y f$ 的加权）

$\begin{aligned} \text{一阶项和} &= 0 + 2 \left( \frac{h}{2} (f_t + f_y f) \right) + 2 \left( \frac{h}{2} (f_t + f_y f) \right) + h (f_t + f_y f) \\ &= h (f_t + f_y f) + h (f_t + f_y f) + h (f_t + f_y f) \\ &= 3h (f_t + f_y f) \end{aligned}$ 代回最终公式中，再乘以外面的 $\frac{h}{6}$： $\frac{h}{6} \left[ 3h (f_t + f_y f) \right] = \mathbf{\frac{h^2}{2} (f_t + f_y f)}$

3. 二阶项（含 $f_{tt}, f_{ty}, f_{yy}$ 和 $f_y(f_t+f_y f)$ 的加权）

这里是整个代数魔术最核心的高潮。我们把 $2k_2$、$2k3$ 和 $k4$ 的二阶项展开相加：

• 对于不含 $f_y$ 的偏导数组合 $\left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} \right]$（我们简记为 $P$）： $\begin{aligned} P \text{ 的系数和} &= 2 \left( \frac{h^2}{8} P \right) + 2 \left( \frac{h^2}{8} P \right) + \frac{h^2}{2} P \\ &= \frac{h^2}{4} P + \frac{h^2}{4} P + \frac{h^2}{2} P \\ &= h^2 P \end{aligned}$
• 对于带有交叉项 $f_y (f_t + f_y f)$ 的部分（我们简记为 $Q$）：
  • $2k_2$ 里面没有这一项（系数为 0）；
  • $2k_3$ 里面含有这一项，系数为 $2 \times \frac{h^2}{8} \times 2 = \frac{h^2}{2}$；
  • $k_4$ 里面含有这一项，系数为 $\frac{h^2}{2}$。 $\begin{aligned} Q \text{ 的系数和} &= 0 + \frac{h^2}{2} Q + \frac{h^2}{2} Q \\ &= h^2 Q \end{aligned}$

把 $P$ 和 $Q$ 合并： $\text{二阶项和} = h^2 (P + Q) = h^2 \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right]$

代回最终公式中，乘以外面的 $\frac{h}{6}$： $\frac{h}{6} \left( h^2 \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right] \right) = \mathbf{\frac{h^3}{6} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right]}$

现在，我们将上面三个部分的加权计算结果全部拼起来，得到 RK4 计算得到的数值解展开： $y_{i+1} = y_i + h f + \frac{h^2}{2} (f_t + f_y f) + \frac{h^3}{6} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right] + O(h^4) \quad \text{--- (RK4 的展开)}$

接着，我们写出真实的微分方程精确解 $y(t_i+h)$ 展开到 $h^3$ 的真实泰勒展开： $y(t_i + h) = y_i + h y'(t_i) + \frac{h^2}{2} y''(t_i) + \frac{h^3}{6} y'''(t_i) + O(h^4)$

我们用全微分（链式法则）求出精确解的各阶导数：

1. 一阶导数： $y'(t_i) = f$
2. 二阶导数： $y''(t_i) = f_t + f_y f$
3. 三阶导数：对 $y'' = f_t + f_y f$ 再次关于时间 $t$ 求导： $\begin{aligned} y'''(t_i) &= \frac{d}{dt} \left[ f_t(t, y) + f_y(t, y) \cdot f(t, y) \right] \\ &= (f_{tt} + f_{ty} y') + (f_{yt} + f_{yy} y') f + f_y (f_t + f_y f) \\ &= f_{tt} + f_{ty} f + (f_{ty} + f_{yy} f) f + f_y (f_t + f_y f) \quad (\text{因 } y'=f, f_{yt}=f_{ty}) \\ &= f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \end{aligned}$

把这些精确导数代入精确解的泰勒展开式： $y(t_i + h) = y_i + h f + \frac{h^2}{2} (f_t + f_y f) + \frac{h^3}{6} \left[ f_{tt} + 2 f f_{ty} + f^2 f_{yy} + f_y (f_t + f_y f) \right] + O(h^4) \quad \text{--- (精确泰勒展开)}$

🏁 惊人的结论：

现在，请将 (RK4 的展开) 和 (精确泰勒展开) 放在一起进行对比。

它们在 $h$ 的零阶项、一阶项、二阶项、三阶项（甚至包含复杂的交叉偏导数 $f_{ty}$ 和 $f_y(f_t+f_y f)$）上，所有的代数项和系数全部一模一样，分毫不差！

这也就意味着，通过这四个斜率的巧妙试探和最终 $1:2:2:1$ 的加权平均：

• 局部截断误差中的 $h^1, h^2, h^3, h^4$ 阶误差全部被完全消去（抵消）。
• 残余的第一项误差出现在 $h^5$ 阶，即局部截断误差为 $O(h^5)$ 。
• 经过区间累积，全局累积误差为 $O(h^4)$。

精度奇迹的代数根源

为什么 RK4 只需要 4 个点就能达到全局四阶精度 ？

• 因为它退化后的辛普森公式本身就是一个具有“代数奇迹”的方法：辛普森公式虽然只用了 $3$ 个点（相当于二次多项式插值），但由于对称性抵消，它对任意三次多项式的积分都是绝对精确的！
• 这种“免费赠送一阶精度”的对称性优势，通过 $k_2$ 和 $k_3$ 在中点的双重采样，被完美遗传给了 RK4。这使得 RK4 拥有极高的代数稳定性，成为了不需要高阶导数就能运行的性价比之王。

局部截断误差： $O(h^5)$，全局误差 $O(h^4)$

🧮 完整数值示例：用 RK4 解 $y' = -2y+1,\ y(0)=1$

import numpy as np
import matplotlib.pyplot as plt

def f(t, y):
    return -2*y + 1

def rk4_step(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h,   y + h * k3)
    return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)

# 参数
y0, h = 1.0, 0.2
t_vals = np.arange(0, 2.01, h)
y_rk4 = [y0]
for i in range(len(t_vals)-1):
    y_rk4.append(rk4_step(f, t_vals[i], y_rk4[-1], h))

# 精确解
y_exact = 0.5 + 0.5 * np.exp(-2 * t_vals)

# 结果对比
print(f"{'t':>6}  {'RK4':>10}  {'Exact':>10}  {'Error':>12}")
for t, yr, ye in zip(t_vals, y_rk4, y_exact):
    print(f"{t:6.2f}  {yr:10.6f}  {ye:10.6f}  {abs(yr-ye):12.2e}")

输出（部分）：

💡 惊人的精度：RK4 用 $h=0.2$，误差约 $5\times 10^{-7}$。Euler 法要达到同样精度需要步长约 $10^{-7}$ 级别！

3.4 自适应步长：Runge-Kutta-Fehlberg (RKF45)

核心想法：每步同时计算 4 阶和 5 阶结果，利用差值估计误差，自动调整步长。

• RKF: 局部截断误差 $O(h^6)$（5阶），可用作误差控制
• Runge-Kutta-Verner: 达到 6 阶

伪代码：

while t < t_end:
    计算 RK4 和 RK5 结果 y4, y5
    误差 err = |y5 - y4|
    if err > tolerance:
        h = h/2; 重算
    else if err < tolerance/10:
        h = 2*h  (放大步长)
    接受结果，前进

💡 物理直觉：解变化慢的地方（平坦区）自动放大步长，变化快的地方（陡峭区）自动缩小步长。就像开车——高速路上踩油门，拐弯处踩刹车。

4. Predictor-Corrector 多步法

4.1 与 Runge-Kutta 的对比

RK4 就是在单步内多点探路 然后代数消去高阶项 达到精度

P-C 预测校正多步法 则是 通过 “不抛弃历史 用插值和积分 重构过去与未来” 来实现 长距离、高效率的稳定计算！！！

RK4 的思路：为了近似这个积分，在当前区间 $[t_i, t_i+h]$ 里临时计算 4 个测试点的斜率 $k_1, k_2, k_3, k_4$，把旧的信息全部扔掉。这需要每前进一步就重复计算 4 次 $f(t,y)$。

多步法的思路：既然我们在前面已经算过了 $t_i, t_{i-1}, t_{i-2}$ 等处的斜率，为什么还要在每一步里重复去算新测试点？我们直接利用这些现成的历史斜率，用多项式把 $f(t, y(t))$ 在过去及当下的曲线“描绘（插值）”出来，然后直接对这个多项式进行精确求积分！

4.2 Adams-Bashforth (显式预测)

三步显式公式：

具体推导

我们要推导三步显式公式： $y_{i+1} = y_i + \frac{h}{12}\left[23f(t_i, y_i) - 16f(t_{i-1}, y_{i-1}) + 5f(t_{i-2}, y_{i-2})\right]$

1. 待定系数法的泰勒展开设定

假设时间步长 $h$ 是恒定的（即 $t_i - t_{i-1} = h$， $t_i - t_{i-2} = 2h$）。为了简洁，我们将已知的历史斜率记为：

• $f_i = f(t_i, y_i)$
• $f_{i-1} = f(t_{i-1}, y_{i-1})$
• $f_{i-2} = f(t_{i-2}, y_{i-2})$

我们希望寻找三个系数 $a, b, c$，使得以下线性组合在代数上具有最高精度： $y(t_{i+1}) \approx y(t_i) + h \left[ a f_i + b f_{i-1} + c f_{i-2} \right] \quad \text{--- (式 1)}$

由于 $y'(t) = f(t, y(t))$，我们把这三项写成导数形式：

• $f_i = y'(t_i)$
• $f_{i-1} = y'(t_{i-1})$
• $f_{i-2} = y'(t_{i-2})$

2. 将所有历史点在当下的 $t_i$ 处进行泰勒展开

我们将 $y(t_{i+1})$，以及历史导数 $y'(t_{i-1})$、$y'(t{i-2})$ 全部在 $ti$ 处展开到 $h^3$ 阶：

1. 左端点 $y(t_{i+1})$ 在 $t_i$ 处的展开（向前走一步 $h$）： $y(t_{i+1}) = y(t_i) + h y'(t_i) + \frac{h^2}{2} y''(t_i) + \frac{h^3}{6} y'''(t_i) + O(h^4) \quad \text{--- (式 2)}$
2. 历史一阶导数 $y'(t_{i-1})$ 在 $t_i$ 处的展开（向后退一步 $-h$）： $y'(t_{i-1}) = y'(t_i) - h y''(t_i) + \frac{h^2}{2} y'''(t_i) + O(h^3) \quad \text{--- (式 3)}$
3. 历史一阶导数 $y'(t_{i-2})$ 在 $t_i$ 处的展开（向后退两步 $-2h$）： $y'(t_{i-2}) = y'(t_i) - (2h) y''(t_i) + \frac{(2h)^2}{2} y'''(t_i) + O(h^3) = y'(t_i) - 2h y''(t_i) + 2h^2 y'''(t_i) + O(h^3) \quad \text{--- (式 4)}$

3. 代入待定系数式，合并同类项

现在，我们将 (式 3) 和 (式 4) 代回我们的待定系数公式 (式 1) 的右端项： $\text{右端项} = y(t_i) + h \left[ a y'(t_i) + b y'(t_{i-1}) + c y'(t_{i-2}) \right]$ 将展开式代入： $\text{右端项} = y(t_i) + h \left[ a y'(t_i) + b \left( y'(t_i) - h y''(t_i) + \frac{h^2}{2} y'''(t_i) \right) + c \left( y'(t_i) - 2h y''(t_i) + 2h^2 y'''(t_i) \right) \right] + O(h^4)$