龙格-库塔

龙格库塔方法（Runge-Kutta Methods） 不是某个人一拍脑袋想出来的，而是经过几十年演化的杰作。让我为你讲述这个精彩的故事和数学推导。

一、历史背景：从欧拉到龙格库塔

1895年：卡尔·龙格（Carl Runge）的突破

当时已知：泰勒展开方法可以得到任意高阶精度 $y_{n + 1} = y_{n} + h y_{n}^{'} + \frac{h ^{2}}{2} y_{n}^{''} + \frac{h ^{3}}{6} y_{n}^{'''} + \dots$
但问题：计算高阶导数 $y^{''}, y^{'''}, \dots$ 非常复杂！
- 对 $y^{'} = f (t, y)$ ，有： $y^{''} = f_{t} + f f_{y}$ $y^{'''} = f_{tt} + 2 f f_{t y} + f^{2} f_{y y} + f_{y} (f_{t} + f f_{y})$ 越来越复杂，尤其是高维情况。

龙格的洞察：能不能不用高阶导数，只用函数值 $f$ 本身，构造高阶方法？

二、核心思想：多”试探”几次

2.1 从几何角度理解

要预测从 $t_{n}$ 到 $t_{n + 1}$ 的路径：

欧拉法：只看起点斜率，直直走过去（容易走偏）
更好的策略：走之前先试探几步，了解”地形”变化

龙格的想法：

在起点算一个斜率 $k_{1}$
用 $k_{1}$ 试探性走到某点
在该点算斜率 $k_{2}$
综合 $k_{1}, k_{2}$ 决定怎么走

三、二阶 RK 方法的推导（龙格，1895年）

3.1 一般形式

设：

y_{n + 1} = y_{n} + h (b_{1} k_{1} + b_{2} k_{2})

其中：

k_{1} = f (t_{n}, y_{n})

k_{2} = f (t_{n} + c_{2} h, y_{n} + a_{21} h k_{1})

待定参数： $b_{1}, b_{2}, c_{2}, a_{21}$ 。

3.2 泰勒展开匹配

精确解展开：

y (t_{n + 1}) = y_{n} + h f + \frac{h ^{2}}{2} (f_{t} + f f_{y}) + O (h^{3})

记 $f = f (t_{n}, y_{n})$ , $f_{t} = \partial f / \partial t$ , $f_{y} = \partial f / \partial y$ 。

数值解展开：
展开 $k_{2}$ （二元泰勒展开）：

k_{2} = f (t_{n} + c_{2} h, y_{n} + a_{21} h k_{1})

= f + c_{2} h f_{t} + a_{21} h f f_{y} + O (h^{2})

代入数值解：

y_{n + 1} = y_{n} + h [b_{1} f + b_{2} (f + c_{2} h f_{t} + a_{21} h f f_{y} + O (h^{2}))]

= y_{n} + h (b_{1} + b_{2}) f + h^{2} b_{2} (c_{2} f_{t} + a_{21} f f_{y}) + O (h^{3})

3.3 匹配系数

与精确解比较：

$h$ 项： $b_{1} + b_{2} = 1$
$h^{2}$ 项： $b_{2} c_{2} = \frac{1}{2}$ , $b_{2} a_{21} = \frac{1}{2}$

三个方程，四个未知数 → 有无穷多解！

3.4 两个著名特解

(1) 中点法（龙格的选择）

取 $c_{2} = \frac{1}{2}$ （走到中点），则：

$b_{2} c_{2} = \frac{1}{2}$ ⇒ $b_{2} = 1$
$b_{1} + b_{2} = 1$ ⇒ $b_{1} = 0$
$b_{2} a_{21} = \frac{1}{2}$ ⇒ $a_{21} = \frac{1}{2}$

公式：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}) y_{n + 1} = y_{n} + h k_{2}

(2) Heun方法（改进欧拉）

取 $c_{2} = 1$ （走到终点），则：

$b_{2} c_{2} = \frac{1}{2}$ ⇒ $b_{2} = \frac{1}{2}$
$b_{1} + b_{2} = 1$ ⇒ $b_{1} = \frac{1}{2}$
$b_{2} a_{21} = \frac{1}{2}$ ⇒ $a_{21} = 1$

公式：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + h, y_{n} + h k_{1}) y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2})

两者都是二阶精度，但系数不同。

四、库塔（Kutta）的贡献（1901年）

4.1 推广到更高阶

威尔海姆·库塔（Wilhelm Kutta）将龙格的思想推广：

用更多中间点（ $k_{3}, k_{4}, \dots$ ）
系统化参数确定方法

一般 s 级 RK 方法：

y_{n + 1} = y_{n} + h i = 1 \sum s b_{i} k_{i}

k_{i} = f (t_{n} + c_{i} h, y_{n} + h j = 1 \sum i - 1 a_{ij} k_{j}), i = 1, \dots, s

参数表（Butcher表）：

c_{1} c_{2} ⋮ c_{s} a_{11} a_{21} ⋮ a_{s 1} b_{1} a_{12} a_{22} ⋮ a_{s 2} b_{2} \dots \dots ⋱ \dots \dots a_{1 s} a_{2 s} ⋮ a_{ss} b_{s}

4.2 经典 RK4 的诞生（库塔，1901年）

库塔系统求解了四阶条件的方程组，得到了最优雅的解：

取对称结构：

$c_{1} = 0, c_{2} = c_{3} = 1/2, c_{4} = 1$
$a_{21} = 1/2, a_{32} = 1/2, a_{43} = 1$ ，其他 $a_{ij} = 0$
$b_{1} = 1/6, b_{2} = 1/3, b_{3} = 1/3, b_{4} = 1/6$

这就是经典的 RK4：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}) k_{3} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}) k_{4} = f (t_{n} + h, y_{n} + h k_{3}) y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

4.3 库塔的数学推导（简化版）

要得到四阶精度，需要匹配泰勒展开到 $h^{4}$ 项。

精确解展开：

y (t_{n + 1}) = y_{n} + h f + \frac{h ^{2}}{2} (f_{t} + f f_{y}) + \frac{h ^{3}}{6} (f_{tt} + 2 f f_{t y} + f^{2} f_{y y} + f_{y} (f_{t} + f f_{y})) + \dots

数值解展开（展开每个 $k_{i}$ 到足够阶，然后组合）会得到很多方程。

库塔巧妙选择对称结构简化方程，然后解出：

从 $\sum b_{i} = 1$
从 $\sum b_{i} c_{i} = 1/2$
从 $\sum b_{i} c_{i}^{2} = 1/3$
从 $\sum b_{i} c_{i}^{3} = 1/4$
以及其他耦合条件…

最终得到经典系数。

五、为什么 RK4 如此优雅？Simpson 法则的联系

5.1 积分视角

将 ODE $y^{'} = f (t, y)$ 写成积分形式：

y (t_{n + 1}) = y (t_{n}) + \int_{t_{n}}^{t_{n + 1}} f (t, y (t)) d t

RK4 本质上是用 Simpson 法则近似这个积分：

Simpson 法则： $\int_{a}^{b} g (t) d t \approx \frac{b - a}{6} [g (a) + 4 g (\frac{a + b}{2}) + g (b)]$
RK4： $\frac{h}{6} [f (起点) + 4 f (中点) + f (终点)]$

但关键：中点的 $f$ 值依赖于解本身！所以需要预测：

用 $k_{1}$ 预测第一个中点值 → $k_{2}$
用 $k_{2}$ 重新预测中点值（更准） → $k_{3}$
用 $k_{3}$ 预测终点值 → $k_{4}$

RK4 = Simpson积分 + 迭代预测修正

六、现代发展：自适应 RK 方法

6.1 嵌入式 RK 方法（Fehlberg, 1969年）

问题：RK4 没有内置误差估计。

解决方案：设计两个不同阶的 RK 方法，共享大部分 $k_{i}$ 计算。

Fehlberg RKF45：

6次函数求值
同时得到 4阶和 5阶结果
用它们的差估计误差

公式（简化）：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + \frac{1}{4} h, y_{n} + \frac{1}{4} h k_{1}) ⋮ k_{6} = f (t_{n} + h, y_{n} + h (- \frac{8}{27} k_{1} + 2 k_{2} - \frac{3544}{2565} k_{3} + \frac{1859}{4104} k_{4} - \frac{11}{40} k_{5}))

4阶结果： $y_{n + 1}^{(4)} = y_{n} + h \sum_{i = 1}^{6} b_{i}^{(4)} k_{i}$
5阶结果： $y_{n + 1}^{(5)} = y_{n} + h \sum_{i = 1}^{6} b_{i}^{(5)} k_{i}$

误差估计： $error \approx ∣ y_{n + 1}^{(5)} - y_{n + 1}^{(4)} ∣$

6.2 Dormand-Prince 方法（1978年）

MATLAB 的 ode45 使用 Dormand-Prince 5(4)：

优化了误差常数
有连续输出插值
成为现代标准

七、RK 方法的数学本质

7.1 级展开（Stage Expansion）

RK 方法本质是用函数值 $f$ 近似高阶导数。

例如，对比泰勒展开：

y_{n + 1} = y_{n} + h f + \frac{h ^{2}}{2} (f_{t} + f f_{y}) + \dots

和 RK2：

y_{n + 1} = y_{n} + h [b_{1} f + b_{2} (f + c_{2} h f_{t} + a_{21} h f f_{y} + \dots)]

通过选择 $b_{2} c_{2} = 1/2, b_{2} a_{21} = 1/2$ ，匹配了 $h^{2}$ 项。

RK 方法：用 $f$ 在不同点的加权组合，模拟泰勒展开中的导数组合。

7.2 树理论（Butcher, 1963年）

约翰·巴特（John Butcher）发展了根树理论，系统化 RK 方法的阶条件。

关键思想：将泰勒展开中的每一项对应到一棵树：

$f$ → 一个节点
$f_{t}$ → 带一个 t-分支的树
$f f_{y}$ → 带一个 y-分支的树
等等

每个树对应一个阶条件方程。

Butcher 表就是这种理论的紧凑表示。

八、为什么 RK 方法如此成功？

8.1 优势

自启动：只需要初值
任意高阶：通过增加级数
无需导数：只计算 $f$
易于实现：公式明确
自适应方便：嵌入式方法

8.2 局限性

计算量大：高阶 RK 需要多次 $f$ 求值
不是 A-稳定：对刚性方程需要极小心
没有记忆效应：不像多步法利用历史信息

九、数值验证：RK4 确实四阶

解 $y^{'} = - y, y (0) = 1$ 到 $t = 1$ ：

$h$	RK4 结果	误差	误差/ $h^{4}$
0.2	0.367914	3.5e-5	0.022
0.1	0.367885	5.5e-6	0.055
0.05	0.367880	8.6e-7	0.055
0.025	0.367879	1.3e-7	0.055

误差 ≈ 常数 × $h^{4}$ ，确认四阶精度。

十、总结：RK 方法的演化路径

时间线：

1895：龙格 → 二阶 RK（避免高阶导数）
1901：库塔 → 系统化，得到 RK4
1963：巴特 → 树理论，完整分类
1969：费尔伯格 → 嵌入式 RK（自适应）
1978：多尔曼德-普林斯 → 优化参数（ode45）

数学路径：

泰勒展开（需要高阶导数）
    ↓ 太复杂！
用多个函数值模拟导数（龙格的洞察）
    ↓
设计权重匹配泰勒系数（库塔的系统化）
    ↓
一般 Butcher 表表示（巴特的理论）
    ↓
嵌入式方法（自适应控制）

核心哲学：用多个点的函数值”询问”微分方程，综合这些信息做出明智的推进决策。

这就是龙格库塔方法：通过聪明的试探与加权，达到了不用计算高阶导数却能获得高阶精度的巧妙平衡。它是数值分析中”用智慧换计算”的典范。

好的！记忆 RK 方法的系数有巧妙的口诀和视觉联想。我来给你一套简单易记的方法。

一、二阶 RK：记住两种常见形式

1. 中点法（Midpoint Method）

公式：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}) y_{n + 1} = y_{n} + h k_{2}

记忆口诀：

“走半步，看斜率，按看到的走全程”

先试探性走半步（ $h /2$ ）
看那里的斜率（ $k_{2}$ ）
按这个斜率走完整的一步（ $h k_{2}$ ）

系数表：

$c_{2} = 1/2$ , $a_{21} = 1/2$
$b_{1} = 0$ , $b_{2} = 1$
简单记：都是 1/2 和 1

2. Heun方法（改进欧拉）

公式：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + h, y_{n} + h k_{1}) y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2})

记忆口诀：

“先猜终点，然后平均”

先用欧拉法猜终点（ $y_{n} + h k_{1}$ ）
算起点和猜的终点的斜率
取平均（ $(k_{1} + k_{2}) /2$ ）

系数表：

$c_{2} = 1$ , $a_{21} = 1$
$b_{1} = 1/2$ , $b_{2} = 1/2$
简单记：除了 b 是 1/2，其他都是 1

二、经典 RK4：记住”1-2-2-1”模式

完整公式：

⎩ ⎨ ⎧ k_{1} = f (t_{n}, y_{n}) k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}) k_{3} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}) k_{4} = f (t_{n} + h, y_{n} + h k_{3}) y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

记忆口诀（四句话）：

一开头，二中前，三中后，四结尾；
一四各一，二三各二，加起来是六。

分解记忆：

1. 四个 k 的计算（时间节点）：

$k_{1}$ ：开头（ $t_{n}$ ）
$k_{2}$ ：中点前（用 $k_{1}$ 预测中点）
$k_{3}$ ：中点后（用 $k_{2}$ 重估中点）
$k_{4}$ ：结尾（用 $k_{3}$ 预测终点）

简单记：1始，2中猜，3中改，4终

2. 四个 k 的参数（步长系数）：

$k_{1}$ ：无（起点）
$k_{2}$ ： $h /2$ （半步）
$k_{3}$ ： $h /2$ （半步）
$k_{4}$ ： $h$ （整步）

简单记：半，半，全（ $k_{1}$ 是 0）

3. 最终权重（关键！）：

\frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

权重模式：1, 2, 2, 1

记忆法：

对称模式：两头小（1），中间大（2）
总和： $1 + 2 + 2 + 1 = 6$ ，所以除以 6

更形象：就像挑担子——中间重（2,2），两头轻（1,1）

视觉记忆法：“梯形+三角形”

画一个区间 ([tn, t{n+1}])，分成两半：

t_n         中点         t_{n+1}
 |------------|------------|
 k1          k2,k3        k4
 1           2,2          1

权重分布：

两个端点（ $k_{1}, k_{4}$ ）：权重 1（像梯形法则的端点）
两个中点（ $k_{2}, k_{3}$ ）：权重 2（合起来相当于中点权重 4，像Simpson法则）

本质：RK4 = Simpson积分法的实现。

快速推导法（考场应急）：

如果忘了系数，可以这样推：

记住 Simpson 法则： $\frac{h}{6} [f_{0} + 4 f_{1/2} + f_{1}]$
但中点值 $f_{1/2}$ 未知，需要预测两次：
- 第一次预测（用 $k_{1}$ ）→ $k_{2}$
- 第二次预测（用 $k_{2}$ ）→ $k_{3}$
把 Simpson 的 4 拆成 2+2： $\frac{h}{6} [k_{1} + 2 k_{2} + 2 k_{3} + k_{4}]$

RK4 弹簧插值

"""
RK4 求解弹簧系统的极简实现
仅使用 Python 标准库，无任何外部依赖
"""
 
import math
from dataclasses import dataclass
from typing import Callable, Optional, List, Tuple
 
 
@dataclass
class SpringState:
    """弹簧状态：位置和速度"""
    x: float  # 位移
    v: float  # 速度
 
 
@dataclass
class SpringDerivative:
    """导数：速度和加速度"""
    dx: float  # dx/dt = v
    dv: float  # dv/dt = a
 
 
class SpringSystem:
    """弹簧-质量-阻尼系统"""
    
    def __init__(self, m: float, k: float, c: float, F: Optional[Callable[[float], float]] = None):
        self.m = m
        self.k = k
        self.c = c
        self.F = F if F else (lambda t: 0.0)
    
    def derivative(self, t: float, y: SpringState) -> SpringDerivative:
        """计算导数 dy/dt"""
        return SpringDerivative(
            dx=y.v,
            dv=(self.F(t) - self.c * y.v - self.k * y.x) / self.m
        )
 
 
def rk4_step(system: SpringSystem, t: float, y: SpringState, h: float) -> SpringState:
    """
    RK4 单步推进
    返回新的状态
    """
    # 四个斜率
    k1 = system.derivative(t, y)
    
    y_mid = SpringState(y.x + h/2 * k1.dx, y.v + h/2 * k1.dv)
    k2 = system.derivative(t + h/2, y_mid)
    
    y_mid = SpringState(y.x + h/2 * k2.dx, y.v + h/2 * k2.dv)
    k3 = system.derivative(t + h/2, y_mid)
    
    y_end = SpringState(y.x + h * k3.dx, y.v + h * k3.dv)
    k4 = system.derivative(t + h, y_end)
    
    # 加权平均更新
    x_new = y.x + h/6 * (k1.dx + 2*k2.dx + 2*k3.dx + k4.dx)
    v_new = y.v + h/6 * (k1.dv + 2*k2.dv + 2*k3.dv + k4.dv)
    
    return SpringState(x_new, v_new)
 
 
def integrate(system: SpringSystem, t0: float, t1: float, y0: SpringState, h: float):
    """
    积分从 t0 到 t1
    返回时间列表和状态列表
    """
    times = [t0]
    states = [y0]
    
    t = t0
    y = y0
    
    n_steps = int((t1 - t0) / h)
    h_actual = (t1 - t0) / n_steps
    
    for _ in range(n_steps):
        y = rk4_step(system, t, y, h_actual)
        t += h_actual
        times.append(t)
        states.append(y)
    
    return times, states

总结：最简记忆

对于 RK4：

四个 k：始、中猜、中改、终
步长：0、半、半、全
权重：1、2、2、1
除以：6（因为 1+2+2+1=6）

对于二阶：

中点法：走半步，看斜率，按它走全程（0,1权重）
Heun法：先猜终点，取平均（1/2,1/2权重）

一句话记 RK4：“始重一中猜二中改三终四，一二二一除六走”

记住这些，你就能在任何需要的时候快速写出 RK 方法的公式了！

Thinking Notes

探索