微分方程计算方法

2026-07-22

📑 目录

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. 微分方程数值解概览

两大问题类型

flowchart TD
    ODE[常微分方程 ODE] --> IVP[初值问题]
    ODE --> BVP[边值问题]
    IVP --> Euler[Euler法 / Taylor展开]
    IVP --> RK[Runge-Kutta族]
    IVP --> PC[Predictor-Corrector]
    BVP --> Shooting[打靶法 Shooting]
    BVP --> FD[有限差分法]
    BVP --> Numerov[Numerov方法]
    Shooting --> Schrodinger[薛定谔方程求解]
    FD --> Linear[转化为线性方程组]

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

2. Euler 法与 Taylor 展开法

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

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

对于一阶 ODE：

$$\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0$$

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

Euler 递推公式：

$$y_{i+1} = y_i + h \cdot f(t_i, y_i)$$

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

📐 数值例子： $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 展开法

递推公式：

$$y_{i+1} = y_i + h y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i + \cdots$$

局部截断误差 (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}]$ 之间 ，多选取

一般形式：

$$ y_{i+1} = y_i + h \sum_{j=1}^{s} b_j k_j $$

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

$$ k_j = f\left(t_i + c_j h,\ y_i + h \sum_{l=1}^{j-1} a_{jl} k_l\right) $$ 这些系数 $a_{jl}$、 $b_j$ 和 $c_j$ 是通过与泰勒级数展开进行复杂的代数比对（Taylor matching）强行解出来

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

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

$$ 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 \\ y_{i+1} = y_i + h \cdot k_2 $$

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

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

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

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

$$ k_1 = f(t_i, y_i) \\ \space \\ k_2 = f(t_i + h,\ y_i + h k_1) \\ \space \\ y_{i+1} = y_i + \frac{h}{2}(k_1 + k_2) $$

说人话：取"起点斜率"和"用 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 的公式里，看看它会退化成什么样：

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$ 居然完全相等了！

$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}$， $f_y = \frac{\partial f}{\partial y}$，且它们都在 $(t_i, 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$

将 上述结果 代入 $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$、$k_2$、$k_3$、$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}$

1. 一阶项（含 $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) \\ \space \\ &= h (f_t + f_y f) + h (f_t + f_y f) + h (f_t + f_y f) \\ \space \\ &= 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)}$

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

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

• 对于不含 $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 \\ \space \\ &= \frac{h^2}{4} P + \frac{h^2}{4} P + \frac{h^2}{2} P \\ \space \\&= 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 \\ \space \\ &= 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]$$

具体推导

我们要推导三步显式公式：

$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})$ 全部在 $t_i$ 处展开到 $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)}$

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) + O(h^3) \quad \text{--- (式 3)}$

1. 历史一阶导数 $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)$$

我们按 $y'(t_i), y''(t_i), y'''(t_i)$ 的阶数进行提取与合并：

$$ \text{右端项} = y(t_i) + h (a + b + c) y'(t_i) + h^2 (-b - 2c) y''(t_i) + h^3 \left( \frac{1}{2}b + 2c \right) y'''(t_i) + O(h^4) \quad \text{--- (式 5)}$$

4. 对比系数，建立线性方程组

为了让数值公式 (式 5) 与真实的泰勒展开 (式 2) 完美契合，我们强制令它们的前三项系数完全相等：

1. $h y'(t_i)$ 的系数： $a + b + c = 1 \quad \text{--- (方程 1)}$
2. $h^2 y''(t_i)$ 的系数： $-b - 2c = \frac{1}{2} \implies b + 2c = -\frac{1}{2} \quad \text{--- (方程 2)}$
3. $h^3 y'''(t_i)$ 的系数： $\frac{1}{2}b + 2c = \frac{1}{6} \quad \text{--- (方程 3)}$

5. 求解方程组，见证系数的诞生

我们用 (方程 2) 减去 (方程 3)： $(b + 2c) - \left(\frac{1}{2}b + 2c\right) = -\frac{1}{2} - \frac{1}{6}$ $\frac{1}{2}b = -\frac{2}{3} \implies b = -\frac{4}{3} = -\frac{16}{12}$

将 $b = -\frac{4}{3}$ 代回 (方程 2)： $-\frac{4}{3} + 2c = -\frac{1}{2} \implies 2c = \frac{4}{3} - \frac{1}{2} = \frac{5}{6} \implies c = \frac{5}{12}$

将 $b = -\frac{16}{12}$ 和 $c = \frac{5}{12}$ 代入 (方程 1)： $a - \frac{16}{12} + \frac{5}{12} = 1 \implies a - \frac{11}{12} = 1 \implies a = \frac{23}{12}$

6. 最终合成

把求得的 $a = \frac{23}{12}, b = -\frac{16}{12}, c = \frac{5}{12}$ 代回 (式 1)，提取公因子 $\frac{h}{12}$：

$y_{i+1} = y_i + \frac{h}{12}\left[23 f_i - 16 f_{i-1} + 5 f_{i-2}\right]$

三步显式 Adams-Bashforth 公式，完美推导完成！

（同理，若将泰勒展开推进到 $O(h^4)$ 阶并引入四个历史点 $f_i, f_{i-1}, f_{i-2}, f_{i-3}$，用完全一样的方法即可推导得出四步显式公式的系数： $\frac{55}{24}, -\frac{59}{24}, \frac{37}{24}, -\frac{9}{24}$。

四步显式公式：

$$y_{i+1} = y_i + \frac{h}{24}\left[55f_i - 59f_{i-1} + 37f_{i-2} - 9f_{i-3}\right]$$

4.3 Adams-Moulton (隐式校正)

隐式：在公式中，包含了 $f(t_{i+1},y_{i+1})$ 这意味着：为了算出下一时刻的状态 $y_{i+1}$ 我们必须在等式右边提前知道下一时刻的状态 $y_{i+1}$ ！！！ 这个在非线性方程中是无法直接进行求解的！

三步隐式公式：

$$y_{i+1} = y_i + \frac{h}{24}\left[9f(t_{i+1}, y_{i+1}) + 19f(t_i, y_i) - 5f(t_{i-1}, y_{i-1}) + f(t_{i-2}, y_{i-2})\right]$$

代数推导

为了推导它的系数，我们假设：

$$ y(t_{i+1}) \approx y(t_i) + h \left[ \alpha f(t_{i+1}, y(t_{i+1})) + \beta f_i + \gamma f_{i-1} + \delta f_{i-2} \right] \quad \text{--- (式 6)} $$ 我们同样把它们写成导数形式，并将所有项在当前点 $t_i$ 处展开到 $h^4$ 阶（因为多引入了一个点 $t_{i+1}$，精度会再提升一阶）：

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) + \frac{h^4}{24} y^{(4)}(t_i) + O(h^5) \quad \text{--- (式 7)}$
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) + \frac{h^3}{6} y^{(4)}(t_i) + O(h^4) \quad \text{--- (式 8)}$
3. 当前一阶导数 $y'(t_i)$： $y'(t_i) = y'(t_i)$
4. 历史一阶导数 $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^{(4)}(t_i) + O(h^4) \quad \text{--- (式 9)}$
5. 历史一阶导数 $y'(t_{i-2})$ 在 $t_i$ 处的展开（向后退 $-2h$）： $y'(t_{i-2}) = y'(t_i) - 2h y''(t_i) + 2h^2 y'''(t_i) - \frac{4h^3}{3} y^{(4)}(t_i) + O(h^4) \quad \text{--- (式 10)}$

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

将 (式 8)、(式 9)、(式 10) 代入 (式 6) 的右侧括号中： $\begin{aligned} \text{右端项} = y(t_i) + h \Big[ &\alpha \left( y'(t_i) + h y''(t_i) + \frac{h^2}{2} y'''(t_i) + \frac{h^3}{6} y^{(4)}(t_i) \right) \\ + &\beta y'(t_i) \\ + &\gamma \left( y'(t_i) - h y''(t_i) + \frac{h^2}{2} y'''(t_i) - \frac{h^3}{6} y^{(4)}(t_i) \right) \\ + &\delta \left( y'(t_i) - 2h y''(t_i) + 2h^2 y'''(t_i) - \frac{4h^3}{3} y^{(4)}(t_i) \right) \Big] \end{aligned}$

按各阶导数项归类合并： $\begin{aligned} \text{右端项} = y(t_i) &+ h (\alpha + \beta + \gamma + \delta) y'(t_i) \\ &+ h^2 (\alpha - \gamma - 2\delta) y''(t_i) \\ &+ h^3 \left( \frac{1}{2}\alpha + \frac{1}{2}\gamma + 2\delta \right) y'''(t_i) \\ &+ h^4 \left( \frac{1}{6}\alpha - \frac{1}{6}\gamma - \frac{4}{3}\delta \right) y^{(4)}(t_i) + O(h^5) \quad \text{--- (式 11)} \end{aligned}$

对比精确展开式 (式 7)，建立方程组

对比 (式 11) 与 (式 7) 的对应系数：

1. $h y'(t_i)$ 对应项： $\alpha + \beta + \gamma + \delta = 1 \quad \text{--- (方程 I)}$
2. $h^2 y''(t_i)$ 对应项： $\alpha - \gamma - 2\delta = \frac{1}{2} \quad \text{--- (方程 II)}$
3. $h^3 y'''(t_i)$ 对应项： $\frac{1}{2}\alpha + \frac{1}{2}\gamma + 2\delta = \frac{1}{6} \implies \alpha + \gamma + 4\delta = \frac{1}{3} \quad \text{--- (方程 III)}$
4. $h^4 y^{(4)}(t_i)$ 对应项： $\frac{1}{6}\alpha - \frac{1}{6}\gamma - \frac{4}{3}\delta = \frac{1}{24} \implies \alpha - \gamma - 8\delta = \frac{1}{4} \quad \text{--- (方程 IV)}$

求解隐式方程组

用 (方程 II) 减去 (方程 IV)： $(\alpha - \gamma - 2\delta) - (\alpha - \gamma - 8\delta) = \frac{1}{2} - \frac{1}{4}$ $6\delta = \frac{1}{4} \implies \delta = \frac{1}{24}$

将 $\delta = \frac{1}{24}$ 代回 (方程 II)： $\alpha - \gamma - \frac{2}{24} = \frac{1}{2} \implies \alpha - \gamma = \frac{14}{24} \quad \text{--- (方程 A)}$

将 $\delta = \frac{1}{24}$ 代入 (方程 III)： $\alpha + \gamma + \frac{4}{24} = \frac{1}{3} \implies \alpha + \gamma = \frac{4}{24} \quad \text{--- (方程 B)}$

联立 (方程 A) 和 (方程 B) 相加： $2\alpha = \frac{18}{24} \implies \alpha = \frac{9}{24}$

联立 (方程 A) 和 (方程 B) 相减： $2\gamma = -\frac{10}{24} \implies \gamma = -\frac{5}{24}$

将已求得的 $\alpha, \gamma, \delta$ 代入 (方程 I) 求 $\beta$： $\frac{9}{24} + \beta - \frac{5}{24} + \frac{1}{24} = 1 \implies \beta + \frac{5}{24} = 1 \implies \beta = \frac{19}{24}$

我们得到了极其完美的系数： $\alpha = \frac{9}{24}, \quad \beta = \frac{19}{24}, \quad \gamma = -\frac{5}{24}, \quad \delta = \frac{1}{24}$

代回 (式 6) 中：

$y_{i+1} = y_i + \frac{h}{24}\left[9f(t_{i+1}, y_{i+1}) + 19f(t_i, y_i) - 5f(t_{i-1}, y_{i-1}) + f(t_{i-2}, y_{i-2})\right]$

三步隐式 Adams-Moulton 公式，完美推导完成！

4.4 PECE 模式（预测-求值-校正-求值）

流程：
1. Predict: 用 Adams-Bashforth (显式) 预测 y_{i+1}^(0)
2. Evaluate: 计算 f(t_{i+1}, y_{i+1}^(0))
3. Correct: 用 Adams-Moulton (隐式) 校正得到 y_{i+1}
4. Evaluate: 重新计算 f(t_{i+1}, y_{i+1})
→ 反复迭代 predictor-corrector 直至收敛

💡 用于 Time-Dependent DFT（含时密度泛函理论）等需要长期稳定积分的高精度计算。

极低的计算开销：RK4 走一步需要计算 4 次极其复杂的 $f(t,y)$。而在 PECE 多步法中，历史斜率都是直接从内存中读取的，走一步通常只需要进行 $1 \sim 2$ 次 $f(t,y)$ 求值（Evaluate 步骤），计算效率比单步法高出 2 到 4 倍。

极强的数值稳定性：隐式方法（Moulton 族）天然具有比显式方法更宽广的绝对稳定域。PECE 通过显式预测和隐式校正的交替，牢牢锁住了非线性长期演化中的能量守恒性，有效防止了长时间积分中的误差发散。

总结

数值分析微分方程中主要的就是两种方法 ：

1、高阶单步法 RK 纵向延伸 不断地k去迭代 只要可以一直算下去 不断精确 消去高阶项

缺点如下

著名的“Butcher障壁”（Butcher's Barrier）：数学家 Butcher 证明了，4 阶 RK 只需要 4步迭代（$k_1 \dots k_4$）；但如果要达到 5 阶精度，最少需要 6 步迭代；要达到 6 阶精度，需要 7 步或 8 步。阶数越高，性价比（精度提升 / 计算量增加）急剧下降。

公式推导极其恐怖：高阶 RK 的系数推导（即代数消去的过程）极其困难，后续往往需要借助计算机代数系统（如 Maple 或 Mathematica）来推导系数。

经典

RK4：性价比之王（4 步达成 4 阶）。

Dormand-Prince (DP54 / solve_ivp 中的 RK45)：通过 7 步迭代（算到 $k_7$），同时拼凑出一个 4 阶和一个 5 阶的公式，用来自动控制步长。

2、P-C多步法 横向时空协作 的 “历史借力”

隐式公式（如 Adams-Moulton）：它的精度极高、稳定性极强（能抗住非线性发散），但它需要解非线性方程（未来点 $y_{i+1}$ 在等式两边都有），非常难算（往往需要牛顿迭代，求导数雅可比矩阵）。

显式公式（如 Adams-Bashforth）：它不含未来点，非常好算，但精度和稳定性稍逊。

==关键==

显式和隐式一结合，显式负责快速给个粗糙的预估，隐式负责精细校正。

最绝的是计算速度：因为在 PECE 模式中，历史的 $f_{i-1}, f_{i-2}$ 早已在之前的步骤中算过并存在内存里了，这步计算它们不需要任何耗时的函数求值（$O(1)$ 读取）。

在长期演化（如你提到的 TDDFT 含时密度泛函、天体轨道计算、分子动力学）中，算一次 $f(t, y)$ 常常要解一次庞大的薛定谔方程或泊松方程。此时，PECE 相比于 RK4，能直接砍掉一半以上的计算时间，同时保持极高的精度！

5. 微分方程组与 Lorenz 混沌

5.1 向量形式

多个变量： $\mathbf{y} = (y_1, y_2, \ldots, y_m)$

方程：

$$\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y})$$

Extended RK4（m 个变量的 RK4）：

for j = 1 to m:   # 每个变量
    k1_j = f_j(t, y_1, y_2, ..., y_m)
    
for j = 1 to m:
    k2_j = f_j(t + h/2, y_1 + h/2·k1_1, ..., y_m + h/2·k1_m)

# 同理计算 k3, k4，然后：
y_j(new) = y_j + h/6 (k1_j + 2k2_j + 2k3_j + k4_j)

在实际写代码时，最容易犯的一个毁灭性错误，就是“一边计算新值，一边覆盖旧值”。

❌ 错误示范：

# 假设我们在算 k2_x, k2_y...

for j in range(m):

    # 危险！当 j=1 (算y) 时，y_mid 需要 x_mid 的预测值。

    # 如果你在这里直接用公式，可能不小心把当前正在被修改的数组当成了输入

    k2[j] = fj [<sup>1</sup>](t + h/2, y[0] + h/2 * k1[0], y[1] + h/2 * k1[1], ...)

正确示范：利用临时的“同步状态向量”

在计算 $k_2, k_3, k_4$ 之前，先用上一个阶段完整的斜率向量，生成一个当前阶段专用的、临时的“探路状态”向量 y_temp，然后所有变量统一在这个 y_temp 下进行函数求值：

Python
# 1. 算完所有变量的 k1

for j in range(m):

    k1[j] = fj [<sup>2</sup>](t, y)

# 2. 同步生成中点状态 1

y_temp = y + (h / 2) * k1  # 这是向量加法，保证所有变量同步向前推进半步

# 3. 在这个同步的中点状态上，算完所有变量的 k2

for j in range(m):

    k2[j] = fj [<sup>3</sup>](t + h / 2, y_temp)

# 4. 同步生成中点状态 2

y_temp = y + (h / 2) * k2

# 5. 算完所有变量的 k3

for j in range(m):

    k3[j] = fj [<sup>3</sup>](t + h / 2, y_temp)

# 6. 同步生成终点状态

y_temp = y + h * k3

# 7. 算完所有变量的 k4

for j in range(m):

    k4[j] = fj [<sup>4</sup>](t + h, y_temp)


# 8. 最终大合流

y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)

💡 注意 $k_2$ 依赖所有变量的 $k_1$——必须算完所有 $k_1$ 再算 $k_2$。

处理多元微分方程组的黄金法则是：“同生共死，步调一致”。 通过引入一个像 y_temp 这样的中间过渡容器，我们在时空上强行将所有互相纠缠的变量锁定在了同一个时刻（起点、中点、终点）。泰勒展开的代数消去魔术也因此得以在多元空间中完美复现，护送整个混沌系统以高精度、高稳定性安全前行。

5.2 Lorenz 方程 — 混沌的经典

$$\begin{aligned} \frac{dx}{dt} &= \sigma(y - x) \\ \frac{dy}{dt} &= x(r - z) - y \\ \frac{dz}{dt} &= xy - bz \end{aligned}$$

参数： $\sigma = 10,\ b = 8/3$，当 $r = 28$ 时出现混沌。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time

# =============================================
# 1. 定义洛伦兹系统（混沌系统）
# =============================================
def lorenz(t, state, sigma=10, r=28, b=8/3):
    """
    洛伦兹方程 (大气对流简化模型)
    dx/dt = sigma*(y - x)
    dy/dt = x*(r - z) - y
    dz/dt = x*y - b*z
    """
    x, y, z = state
    dx = sigma * (y - x)
    dy = x * (r - z) - y
    dz = x * y - b * z
    return np.array([dx, dy, dz])

# =============================================
# 2. 核心：向量版 RK4 步进函数（通用，支持任意维度的系统）
# =============================================
def rk4_step_vector(f, t, y, h):
    """
    四阶龙格-库塔法 (向量形式)
    参数:
        f: 微分方程函数，返回 np.array
        t: 当前时间 (标量)
        y: 当前状态向量 (np.array)
        h: 步长
    返回:
        下一个时间步的状态向量 (np.array)
    """
    # 四个斜率计算（全为向量运算）
    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)
    
    # 加权平均更新状态
    y_next = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
    return y_next

# =============================================
# 3. 主程序：求解并保存轨迹
# =============================================
if __name__ == "__main__":
    # 参数设置
    state = np.array([1.0, 1.0, 1.0])  # 初始条件 (x0, y0, z0)
    h = 0.01                           # 步长 (越小精度越高，但计算量越大)
    N = 5000                           # 迭代步数 (总积分时间 T = h*N = 50秒)
    
    # 记录轨迹（为了节省内存，可以预先分配数组，但用列表 append 更直观）
    trajectory = [state.copy()]
    t_current = 0.0
    
    print(f"开始求解洛伦兹系统，步长 h={h}，总步数 N={N}...")
    start_time = time.time()
    
    # 主循环
    for i in range(N):
        state = rk4_step_vector(lorenz, t_current, state, h)
        t_current += h
        trajectory.append(state.copy())  # 注意要 copy，防止后续修改影响历史数据
        
        # 可选：每 500 步打印一次进度
        if (i + 1) % 500 == 0:
            print(f"  已完成 {i+1}/{N} 步，当前时间 t = {t_current:.2f}")
    
    end_time = time.time()
    print(f"计算完成！耗时 {end_time - start_time:.4f} 秒。")
    
    # 将轨迹列表转换为 numpy 数组，方便切片绘图
    traj = np.array(trajectory)
    print(f"轨迹数组形状: {traj.shape}")  # 应为 (N+1, 3)

    # =============================================
    # 4. 可视化：绘制三维相图和二维投影
    # =============================================
    fig = plt.figure(figsize=(15, 5))
    
    # 子图 1: 三维相空间 (经典蝴蝶效应)
    ax1 = fig.add_subplot(131, projection='3d')
    ax1.plot(traj[:, 0], traj[:, 1], traj[:, 2], linewidth=0.8, color='blue')
    ax1.set_xlabel('X')
    ax1.set_ylabel('Y')
    ax1.set_zlabel('Z')
    ax1.set_title('洛伦兹吸引子 (3D)')
    
    # 子图 2: X-Z 投影 (最常见的蝴蝶图)
    ax2 = fig.add_subplot(132)
    ax2.plot(traj[:, 0], traj[:, 2], linewidth=0.6, color='red')
    ax2.set_xlabel('X')
    ax2.set_ylabel('Z')
    ax2.set_title('X-Z 投影')
    ax2.grid(True, linestyle='--', alpha=0.5)
    
    # 子图 3: X-Y 投影
    ax3 = fig.add_subplot(133)
    ax3.plot(traj[:, 0], traj[:, 1], linewidth=0.6, color='green')
    ax3.set_xlabel('X')
    ax3.set_ylabel('Y')
    ax3.set_title('X-Y 投影')
    ax3.grid(True, linestyle='--', alpha=0.5)
    
    plt.tight_layout()
    plt.show()
    
    # =============================================
    # 5. 额外：打印最终状态和简单的统计学检查
    # =============================================
    final_state = traj[-1]
    print(f"\n最终状态 (t={t_current:.2f}): X={final_state[0]:.6f}, Y={final_state[1]:.6f}, Z={final_state[2]:.6f}")
    
    # 检查是否有数值溢出 (洛伦兹系统在参数 28 时是混沌有界的，一般不会溢出)
    if np.any(np.isnan(traj)):
        print("⚠️ 警告：检测到 NaN，请尝试减小步长 h 或检查初始条件！")
    else:
        print("✅ 数值稳定，未出现溢出。")

💡 混沌标志：初始条件的微小差异（如 $x_0 = 1.000$ vs $1.001$）在 $t \approx 15$ 后演化出完全不同的轨迹——蝴蝶效应！

6. 高阶微分方程的降阶处理

6.1 核心技巧：高阶 → 一阶方程组

二阶 ODE： $y'' + p(t)y' + q(t)y = r(t)$

定义新变量降阶：

$$\begin{aligned} u_1 &= y \\ u_2 &= y' \end{aligned} \quad\Rightarrow\quad \begin{cases} u_1' = u_2 \\ u_2' = r(t) - p(t)u_2 - q(t)u_1 \end{cases}$$

这就是 $m$ 个一阶方程 → 用 Extended RK4 直接求解。

6.2 📐 完整例子：阻尼谐振子

方程： $y'' + 2\gamma y' + \omega_0^2 y = 0$，设 $\gamma=0.3,\ \omega_0=2$，$y(0)=1,\ y'(0)=0$

降阶为：

$$\begin{cases} u_1' = u_2 \\ u_2' = -0.6u_2 - 4u_1 \end{cases}$$

def damped_oscillator(t, u):
    u1, u2 = u
    gamma, omega0 = 0.3, 2.0
    du1 = u2
    du2 = -2*gamma*u2 - omega0**2 * u1
    return np.array([du1, du2])

u0 = np.array([1.0, 0.0])  # [y(0), y'(0)]
h = 0.05
t_vals = np.arange(0, 10.0, h)
u_traj = [u0]
for i in range(len(t_vals)-1):
    u_traj.append(rk4_step_vector(damped_oscillator, t_vals[i], u_traj[-1], h))

# 结果：振幅随时间指数衰减，频率略小于 ω₀

💡 通用策略：任意 $n$ 阶 ODE 都可以转化为 $n$ 个一阶方程。这个技巧在计算物理中无处不在。

7. 刚性微分方程（Stiff ODE）

7.1 什么是刚性？

系统中同时存在"快变量"和"慢变量"时，显式方法必须用极小的步长来追踪最快分量，否则发散。

快变量：变化速度极快 在瞬间 衰减到0

慢变量： 真正关心的 主导了系统长期演变的变量

刚性的本质：微分方程组的 Jacobi 矩阵特征值 $\lambda_i$ 中，$\max|\lambda_i| / \min|\lambda_i|$ 很大。

7.2 稳定性对比

致命的痛点：为什么显式方法（如欧拉、RK4）会在这里死掉？

如果你想计算这个系统在 $t = 5$ 秒时的状态。此时弹簧 A 早就静止了，照理说我们应该可以跨大步子（比如步长 $h = 0.1$ 秒）来快速计算弹簧 B 的悠闲运动。

然而，如果你使用显式欧拉法（或者 RK4）： 只要你的步长 $h$ 稍稍大于 0.0001 秒（比如你用了 $h = 0.01$），那个早就已经静止的“快变量”就会在数值计算中像疯了一样成倍放大、瞬间发散到无穷大！

显式方法被逼着必须使用 $h = 0.00001$ 秒的超级小步长，去算那毫无意义、早就静止了的快变量，从而在极小的步长上把自己硬生生“耗死”。

绝对稳定性与特征值：数值世界里的“紧箍咒”

为了看清为什么显式方法会发散，而隐式方法不会，我们必须做一次严密的代数推导。

1. 显式欧拉的自我毁灭推导

假设我们面对最简单的一个一阶线性微分方程（代表我们系统里的“快变量”衰减过程）： $\frac{dy}{dt} = \lambda y, \quad (\text{其中 } \lambda < 0 \text{，是一个极大的负数，比如 } \lambda = -1000)$ 它的精确解析解是 $y(t) = C e^{\lambda t} = C e^{-1000 t}$。随着时间推移，它会以极快的速度衰减到 $0$。

现在，我们用显式欧拉法来计算它： $y_{i+1} = y_i + h f(t_i, y_i) = y_i + h (\lambda y_i)$ 提取公因子 $y_i$： $y_{i+1} = (1 + h\lambda) y_i \quad \text{--- (式 1)}$

致命的死亡条件：

要想让计算出来的 $y$ 像真实物理世界一样不断衰减（不发散），我们必须保证每前进一步，新值都比旧值小。也就是： $|1 + h\lambda| < 1$ 由于 $\lambda = -1000$，我们解这个绝对值不等式： $-1 < 1 - 1000h < 1 \implies 0 < 1000h < 2 \implies h < 0.002$

瞧见了吗？不管你有多想跨大步子，你的步长 $h$ 都被死死地限制在了 $0.002$ 秒以内！ 一旦你的 $h$ 哪怕大了一丁点（比如 $h = 0.003$），系数就会变成 $|1 + h\lambda| = |1 - 3| = 2$。 下一步的值就会变成上一步的 2 倍，再下一步是 4 倍、8 倍……瞬间发散崩溃！

这就是表里的： $\text{Explicit: } |1 + h\lambda| < 1$。

隐式梯形法的“降维打击”

现在，我们把相同的方程，用隐式梯形法来算。 隐式梯形法的公式为： $y_{i+1} = y_i + \frac{h}{2} [f(t_i, y_i) + f(t_{i+1}, y_{i+1})]$

代入 $f(t,y) = \lambda y$： $y_{i+1} = y_i + \frac{h}{2} [\lambda y_i + \lambda y_{i+1}]$

我们将含有未来项 $y_{i+1}$ 的项移到等式的左边： $y_{i+1} - \frac{h\lambda}{2} y_{i+1} = y_i + \frac{h\lambda}{2} y_i$ 提取并整理： $\left(1 - \frac{h\lambda}{2}\right) y_{i+1} = \left(1 + \frac{h\lambda}{2}\right) y_i$

$y_{i+1} = \left( \frac{1 + \frac{h\lambda}{2}}{1 - \frac{h\lambda}{2}} \right) y_i \quad \text{--- (式 2)}$

🎉 奇迹发生了：

我们要让数值解不发散，同样需要放大系数的绝对值小于 1： $\left| \frac{1 + \frac{h\lambda}{2}}{1 - \frac{h\lambda}{2}} \right| < 1$

由于物理衰减中 $\lambda < 0$（即负实部），所以：

• 分母 $1 - \frac{h\lambda}{2}$ 一定是一个大于 1 的正数。
• 分子 $1 + \frac{h\lambda}{2}$ 的绝对值，在 $\lambda < 0$ 时，永远小于分母！

这意味着：无论你的步长 $h$ 取得有多大（哪怕 $h = 10000$ 甚至无穷大），这个放大系数的绝对值都永远小于 1！ 隐式梯形法是无条件稳定的。你可以放心大胆地用巨无霸步长跨过去，那个快变量只会温顺地、迅速地衰减到 0，绝对不会发散。

7.3 隐式梯形法 (Implicit Trapezoidal)

$$y_{i+1} = y_i + \frac{h}{2}\left[f(t_i, y_i) + f(t_{i+1}, y_{i+1})\right]$$

隐式方法虽然好，但是它有一个致命难点：

未来项 $y_{i+1}$ 藏在复杂的非线性函数 $f(t_{i+1}, y_{i+1})$ 里面，我们该怎么把它解出来？

以隐式梯形法为例，我们定义一个残差函数 $F(y_{i+1})$，我们的目标是找到一个 $y_{i+1}$ 使得这个函数等于 0： $F(y_{i+1}) = y_{i+1} - y_i - \frac{h}{2} \left[ f(t_i, y_i) + f(t_{i+1}, y_{i+1}) \right] = 0$

因为 $f$ 是非线性的，我们必须在每一步里，都用牛顿迭代法（Newton-Raphson）去把它活生生“猜”出来。

牛顿迭代手算推导：

牛顿迭代法是求解非线性方程 $F(x) = 0$ 最犀利的武器。它的几何直觉极其美妙：

1. 第一步（盲猜）：我们先闭着眼睛对 $y_{i+1}$ 猜一个初始值，记为 $y^{(0)}$。

• 注：怎么猜？最常用的办法是用最简单的显式欧拉法跨一步： $y^{(0)} = y_i + h f(t_i, y_i)$，这个粗糙的值作为我们的起点。

1. 第二步（作切线）：我们画出函数 $F(x)$ 的曲线，找到曲线在 $x = y^{(0)}$ 处的点 $(y^{(0)}, F(y^{(0)}))$。我们在这个点上作一条切线。
2. 第三步（切线求交点）：这条切线的斜率就是导数 $F'(y^{(0)})$。切线与 $x$ 轴（即 $F(x) = 0$ 处）相交的点，就会比我们盲猜的 $y^{(0)}$ 更接近真正的根。我们将这个交点记为下一次的猜测值 $y^{(1)}$。

我们用一次函数（切线方程）来写出这个过程： 在点 $x = y^{(k)}$ 处的切线方程为： $y - F(y^{(k)}) = F'(y^{(k)}) (x - y^{(k)})$ 我们要找它与 $x$ 轴的交点，即令 $y = 0$，求解此时的 $x$（即下一次逼近值 $y^{(k+1)}$：

$（0 - F(y^{(k)}) = F'(y^{(k)}) (y^{(k+1)} - y^{(k)})$

移项整理： $y^{(k+1)} - y^{(k)} = -\frac{F(y^{(k)})}{F'(y^{(k)})}$

$y^{(k+1)} = y^{(k)} - \frac{F(y^{(k)})}{F'(y^{(k)})} \quad \text{--- (牛顿迭代核心公式)}$

我们要对 $F(y_{i+1})$ 关于 $y_{i+1}$ 求导： $F'(y_{i+1}) = \frac{\partial}{\partial y_{i+1}} \left( y_{i+1} - y_i - \frac{h}{2} [f(t_i, y_i) + f(t_{i+1}, y_{i+1})] \right) = 1 - \frac{h}{2} \frac{\partial f}{\partial y}$ 这里的 $\frac{\partial f}{\partial y}$ 就是一维时的导数，在多元微分方程组里，它就是赫赫有名的雅可比矩阵（Jacobian Matrix） $J$。

因此，牛顿迭代的每一次修正公式为： $y^{(k+1)} = y^{(k)} - \frac{y^{(k)} - y_i - \frac{h}{2} [f(t_i, y_i) + f(t_{i+1}, y^{(k)})]}{1 - \frac{h}{2} J}$

$y_{i+1}$ 出现在方程右侧，需用 Newton 迭代求解：

# 伪代码：隐式梯形法一步
def trapezoidal_step(f, t, y, h):
    # 先用 Euler 预测一个初猜
    y_guess = y + h * f(t, y)
    
    # Newton 迭代求解
    for _ in range(max_iter):
        F = y_guess - y - h/2 * (f(t, y) + f(t+h, y_guess))
        if abs(F) < tol: break
        # J = 1 - h/2 * df/dy (对单变量)
        J = 1 - h/2 * jacobian(f, t+h, y_guess)
        y_guess = y_guess - F/J
    
    return y_guess

7.4 Gear 方法 (2阶 / 3阶)

Backward Differentiation Formula

普通多步法（Adams 族）：通过对斜率 $f$ 进行插值，去近似积分。

Gear 方法（BDF 族）：通过对状态 $y$（即当前点和历史点）进行高阶多项式插值，强行对这个多项式求导，以此来逼近未来的斜率。

2 阶 Gear：

$$y_{i+1} = \frac{4}{3}y_i - \frac{1}{3}y_{i-1} + \frac{2h}{3}f(t_{i+1}, y_{i+1})$$

它的代数原理极其简单： 假设我们在 $t_{i+1}, t_i, t_{i-1}$ 这三个点上，穿过状态值 $y_{i+1}, y_i, y_{i-1}$ 作一条二次抛物线 $P(t)$。 通过拉格朗日插值法，我们可以写出这条代表位移的抛物线 $P(t)$ 的表达式。 然后，我们强制要求抛物线在未来时刻 $t_{i+1}$ 处的导数 $P'(t_{i+1})$，必须完美等于物理方程给出的斜率 $f(t_{i+1}, y_{i+1})$！ 即： $P'(t_{i+1}) = f(t_{i+1}, y_{i+1})$ 对插值多项式求导并整理后，这一步就会硬生生、分毫不差地凑出： $y_{i+1} = \frac{4}{3}y_i - \frac{1}{3}y_{i-1} + \frac{2h}{3} f(t_{i+1}, y_{i+1})$

同理，若用 4 个点作三次抛物线，求导后就能推导出 3 阶 Gear：

$y_{i+1} = \frac{18}{11}y_i - \frac{9}{11}y_{i-1} + \frac{2}{11}y_{i-2} + \frac{6h}{11}f(t_{i+1}, y_{i+1})$

3 阶 Gear：

$$y_{i+1} = \frac{18}{11}y_i - \frac{9}{11}y_{i-1} + \frac{2}{11}y_{i-2} + \frac{6h}{11}f(t_{i+1}, y_{i+1})$$

都是隐式方法，每步需 Newton 迭代求解 $y_{i+1}$。

💡 稳定性准则（微分方程组）：

$$\text{Explicit: } |1 + h\lambda_i| < 1 \quad \text{vs} \quad \text{Implicit: 无条件稳定（$Re(\lambda_i)<0$ 时）}$$

为什么 Gear 方法（BDF）最适合 Stiff 方程？（深入绝对稳定性）

为了彻底看清为什么这套方法能成为刚性（Stiff）微分方程的“降维打击”王牌，我们需要从绝对稳定域的数学本源去剖析它。

1. 测试方程与线性化

为了分析任何一个数值方法的稳定性，数学家都会使用标准线性测试方程： $$ \frac{dy}{dt} = \lambda y \quad (\text{其中 } \lambda \in \mathbb{C}, \text{ 复数，且其实部 } \text{Re}(\lambda) < 0) $$ 这个方程代表了物理系统中的衰减分量（快变量）。我们定义无量纲参数 $z = h\lambda$。 在真实物理世界里，因为实部为负，所以随着 $t \to \infty$，解析解 $y(t) \to 0$。 我们要求：数值方法计算出来的序列 $y_n$ 在步数 $n \to \infty$ 时，也必须衰减到 $0$。

2. BDF2 的特征方程与绝对稳定条件

我们将测试方程 $f(t, y) = \lambda y$ 代入刚刚推导出的 BDF2 公式中： $$ y_{i+1} = \frac{4}{3} y_i - \frac{1}{3} y_{i-1} + \frac{2}{3} z y_{i+1} \quad (\text{其中 } z = h\lambda) $$ 我们将所有项移到左侧并整理： $$ \left( 1 - \frac{2}{3} z \right) y_{i+1} - \frac{4}{3} y_i + \frac{1}{3} y_{i-1} = 0 $$ 为了去分母，两边乘以 $3$： $$ (3 - 2z) y_{i+1} - 4 y_i + y_{i-1} = 0 $$ 这是一个关于序列 $y$ 的二阶常系数线性差分方程。其通解的形式为 $y_n = C r^n$，其中 $r$ 是其特征方程的根。 我们将 $y_n = r^n$ 代入上式，得出其特征方程： $$ (3 - 2z) r^2 - 4 r + 1 = 0 $$ 根据线性差分方程的稳定性理论，要让数值解 $y_n \to 0$ 稳定不发散，特征方程的两个根 $r_1, r_2$ 的模长都必须严格小于 1（即它们必须落在复平面的单位圆内）： $$ |r_1| < 1 \quad \text{且} \quad |r_2| < 1 $$ 我们来解出特征根 $r$ 看看它与 $z = h\lambda$ 的关系： $$ r = \frac{4 \pm \sqrt{16 - 4(3 - 2z)}}{2(3 - 2z)} = \frac{4 \pm \sqrt{4 + 8z}}{2(3 - 2z)} = \frac{2 \pm \sqrt{1 + 2z}}{3 - 2z} $$

🌟 为什么它极其稳定？

假设我们的方程存在极度刚性的快分量，这意味着 $\lambda$ 是一个极大的负实数（例如 $\lambda = -10000$）。如果我们用大步长 $h = 1$ 去跨，那么 $z = h\lambda = -10000$ 就会变成一个极其恐怖的负数。 我们来看一下当 $z \to -\infty$ 时，特征根 $r$ 会发生什么： $$ \lim_{z \to -\infty} r = \lim_{z \to -\infty} \frac{2 \pm \sqrt{1 + 2z}}{3 - 2z} $$ 由于分母上有 $z$ 的一次项，而分子上只有根号下 $z$（即半次方项）。根据极限法则，分母的增长速度远快于分子： $$ \lim_{z \to -\infty} r = 0 $$ 这意味着什么？

当面对极强刚性（ $z \to -\infty$）时，BDF2 的特征根 $r$ 直接趋近于 0！ 特征根趋近于 0，意味着每走一步，误差不仅不会放大，反而会被以近乎 100% 的效率瞬间抹除、强行掐灭！ 这就是为什么说它是“黑洞一样的结构，能够吞噬快分量误差”的数学解释。

3. 绝对稳定域的对比（显式 vs 隐式）

为了让你更具象地感受，我们画出复平面上的“绝对稳定域”（也就是能让特征根模长小于 1 的所有 $z$ 的集合）：

• 显式欧拉法：稳定域是复平面上以 $(-1, 0)$ 为圆心、半径为 $1$ 的小圆圈。 $z = h\lambda$ 必须在这个小圆圈内。当 $\lambda = -1000$ 时， $h$ 必须小于 $0.002$ 才能把 $z$ 关进这个圆圈。
• 隐式梯形法：稳定域是整个复平面的左半平面（A-Stable）。只要 $\text{Re}(\lambda) < 0$，无论 $h$ 有多大， $z$ 都永远在稳定域里。但是当 $z \to -\infty$ 时，其放大系数趋于 $-1$（边界上），误差虽然不发散，但会产生持久的数值振荡。
• Gear 方法（BDF2）：不仅是 A-stable 的（整个左半平面几乎全包），而且它拥有极强的 L-稳定性（L-Stable）。也就是说，当 $z \to -\infty$ 时，它的放大系数直接收敛到 0。这相当于对快分量的数值误差实施了“降维打击”——只要你刚性发作得越厉害，它掐死误差的速度就越快。

核心结论：Implicit 的稳定性远优于 Explicit。对 Stiff 问题，宁可用隐式方法多算几步 Newton 迭代，也别用显式方法在极小步长上耗死自己。

为什么 Gear 方法最适合 Stiff 方程？

因为它的左边是纯粹的状态项 $y$，右边只有一项未来的斜率项 $f(t_{i+1}, y_{i+1})$。 在代数合并和牛顿迭代时，这种结构在复平面上的绝对稳定区域极其庞大，能够像黑洞一样，把所有具有极大负实部特征值（也就是死得最快的那些快分量）的误差全部吞噬、彻底压制。

8. 边值问题BVP与打靶法（Shooting）

8.1 问题描述

边值问题：已知 $y(a) = \alpha,\ y(b) = \beta$，求解 $y(x)$。

与初值问题的根本不同——两端都有条件。

核心打靶法就是：

把未知的初始斜率 $$ t = y'(a) $$ 当成一个可调参数，不断试，直到积分到 $b$ 时满足 $$ y(b)=\beta $$ 打靶法核心思路：把边值问题转化为初值问题！

8.2 线性打靶法

优势之处：线性方程满足叠加原理

对线性方程： $$ y'' + p(x)y' + q(x)y = r(x), \quad y(a)=\alpha,\ y(b)=\beta $$ 步骤：

1. 解两个辅助初值问题：

$$\begin{aligned} \text{问题1: } & y_1'' + py_1' + qy_1 = r, \quad y_1(a)=\alpha,\ y_1'(a)=0 \\ \text{问题2: } & y_2'' + py_2' + qy_2 = 0, \quad y_2(a)=0,\ y_2'(a)=1 \end{aligned}$$

1. 用 RK4 分别求 $y_1(b),\ y_2(b)$

1. 混合两个解： $y(x) = y_1(x) + t \cdot y_2(x)$

其中系数 $t$ 由边界条件确定：

$$t = \frac{\beta - y_1(b)}{y_2(b)}$$

两个解都满足方程（因为是线性叠加），调整混合比例 $t$ 来同时满足两端边界条件。

线性打靶法本质

你可以理解为： $$ y_1 $$ 是“先打一枪，斜率取 0”的轨迹。 $$ y_2 $$ 是“单位斜率造成的响应”。

然后真正需要的斜率修正量就是： $$ t=\frac{\text{终点缺口}}{\text{单位斜率带来的终点变化}} $$ 即： $$ t = \frac{\beta-y_1(b)}{y_2(b)} $$

注意：什么时候会失败？

如果 $$ y_2(b)=0 $$ 那么公式分母为零。

这说明初始斜率的变化对终点 $y(b)$ 没有影响，或者对应的齐次边值问题存在非零解。

这通常意味着：

• 问题无解；
• 或者解不唯一；
• 或者参数正好落在某个本征值上。

例如 Sturm-Liouville 本征值问题中，齐次解可以同时满足两端边界，这时 Green 函数分母 Wronskian 也会出问题。

8.3 非线性打靶法 + 割线法求 $t$

对于非线性问题： $$ y''=f(x,y,y') $$ 边界条件： $$ y(a)=\alpha,\qquad y(b)=\beta $$ 我们仍然设： $$ t=y'(a) $$ 每给一个 $t$，从 $a$ 积分到 $b$，得到： $$ y(b;t) $$ 定义残差函数： $$ F(t)=y(b;t)-\beta $$ 目标： $$ F(t)=0 $$

割线法推导

假设已经有两个试探斜率： $$ t_{k-1},\qquad t_k $$ 对应残差：

$$ F(t_{k-1}),\qquad F(t_k) $$ 我们用一条直线近似 $F(t)$。

割线斜率是： $$ \frac{F(t_k)-F(t_{k-1})}{t_k-t_{k-1}} $$ 直线方程为： $$ F(t)\approx F(t_k) + \frac{F(t_k)-F(t_{k-1})}{t_k-t_{k-1}}(t-t_k) $$ 我们希望下一次 $t_{k+1}$ 使得这个线性近似为 0： $$ 0= F(t_k) + \frac{F(t_k)-F(t_{k-1})}{t_k-t_{k-1}}(t_{k+1}-t_k) $$ 移项： $$ \frac{F(t_k)-F(t_{k-1})}{t_k-t_{k-1}}(t_{k+1}-t_k) = -F(t_k) $$ 所以 $$ t_{k+1}-t_k = -F(t_k) \frac{t_k-t_{k-1}}{F(t_k)-F(t_{k-1})} $$ 得到： $$ t_{k+1} = t_k - F(t_k) \frac{t_k-t_{k-1}}{F(t_k)-F(t_{k-1})} $$ 由于 $$ F(t)=y(b;t)-\beta $$ 所以 $$ t_{k+1} = t_k - \bigl(y(b;t_k)-\beta\bigr) \frac{t_k-t_{k-1}} {y(b;t_k)-y(b;t_{k-1})} $$ 这就是你写的公式。

非线性打靶法算法

1. 选两个初始猜测：

$$ t_0,\qquad t_1 $$

1. 用 RK4 积分：

$$ y(b;t_0),\qquad y(b;t_1) $$

1. 计算残差：

$$ F(t_0)=y(b;t_0)-\beta $$

1. 用割线法更新：

$$ t_2= t_1 - F(t_1) \frac{t_1-t_0}{F(t_1)-F(t_0)} $$

1. 重复直到：

$$ |F(t_k)|=|y(b;t_k)-\beta|<\varepsilon $$

💡 为什么叫"打靶"：你站在 a 点，调整出射角 $t = y'(a)$，"射击"目标是 b 点的边界值。射偏了就调整角度再来一枪，直到命中靶心。

8.4 割线法与牛顿法的选择

如果能算出： $$ F'(t)=\frac{\partial y(b;t)}{\partial t} $$ 那么可以用牛顿法： $$ t_{k+1} = t_k-\frac{F(t_k)}{F'(t_k)} $$

对打靶法，$dy(b,t)/dt$ 往往不好求，割线法更实用。

9. Numerov 方法

9.1 专门为二阶 ODE 设计

Numerov 方法针对以下形式的方程（无 $y'$ 项）： $$y'' + k(x)y = F(x)$$

普通二阶方程： $$ y''+p(x)y'+q(x)y=r(x) $$ 包含 $y'$，Numerov 不能直接套。

但有时可以通过变量替换消去一阶导。

令： $$ y(x)=u(x)e^{-\frac12\int p(x)\,dx} $$ 设 $$ S(x)=\frac12\int p(x)\,dx $$ 所以： $$ y=e^{-S}u $$ 求导： $$ y'=e^{-S}(u'-S'u) $$ 再求导： $$ y'' = e^{-S} \left[ u''-2S'u'+(S'^2-S'')u \right] $$ 代入原方程： $$ y''+p y'+q y=r $$ 得到： $$ e^{-S} \left[ u''-2S'u'+(S'^2-S'')u \right] + p e^{-S}(u'-S'u) + q e^{-S}u=r $$ 两边乘以 $e^S$： $$ u''-2S'u'+(S'^2-S'')u + p u'-pS'u+q u = e^S r $$ 因为 $$ S'=\frac p2 $$ 所以一阶导项系数： $$ -2S'+p=-p+p=0 $$ 一阶导消失。

剩下 $u$ 的系数： $$ S'^2-S''-pS'+q $$ 代入： $$ S'^2=\frac{p^2}{4} $$ 所以： $$ S'^2-S''-pS'+q = \frac{p^2}{4}-\frac{p'}{2}-\frac{p^2}{2}+q $$ 因此变成： $$ u'' + \left( q-\frac{p'}{2}-\frac{p^2}{4} \right)u = e^S r $$ 这就变成了 Numerov 形式： $$ u''+K(x)u=\tilde F(x) $$ 其中： $$ K(x)=q-\frac{p'}{2}-\frac{p^2}{4} $$

9.2 递推公式

我们从 Taylor 展开开始。

在网格： $$ x_i=x_0+ih $$ 记： $$ y_i=y(x_i),\qquad k_i=k(x_i),\qquad F_i=F(x_i) $$ 方程： $$ y''+ky=F $$ 也就是： $$ y''=F-ky $$

第一步：中心差分展开

Taylor 展开： $$ y_{i+1} = y_i+h y_i' +\frac{h^2}{2}y_i'' +\frac{h^3}{6}y_i''' +\frac{h^4}{24}y_i^{(4)} +\frac{h^5}{120}y_i^{(5)} +\frac{h^6}{720}y_i^{(6)} +\cdots $$ 相加： $$ y_{i+1}+y_{i-1} = 2y_i +h^2 y_i'' +\frac{h^4}{12}y_i^{(4)} +\frac{h^6}{360}y_i^{(6)} +\cdots $$ 所以： $$ y_{i+1}-2y_i+y_{i-1} = h^2y_i'' +\frac{h^4}{12}y_i^{(4)} +O(h^6) $$

第二步：用差分近似 $y^{(4)}$

因为： $$ y^{(4)}=(y'')'' $$ 所以可以用二阶中心差分近似： $$ y_i^{(4)} \approx \frac{y_{i+1}''-2y_i''+y_{i-1}''}{h^2} $$ 代入： $$ y_{i+1}-2y_i+y_{i-1} = h^2 y_i'' + \frac{h^4}{12} \frac{y_{i+1}''-2y_i''+y_{i-1}''}{h^2} +O(h^6) $$ 整理： $$ y_{i+1}-2y_i+y_{i-1} = h^2 \left[ y_i'' + \frac{1}{12} \left( y_{i+1}''-2y_i''+y_{i-1}'' \right) \right] +O(h^6) $$ 把括号内合并： $$ y_i'' + \frac{1}{12} (y_{i+1}''-2y_i''+y_{i-1}'') $$ 所以： $$ y_{i+1}-2y_i+y_{i-1} = \frac{h^2}{12} \left( y_{i+1}''+10y_i''+y_{i-1}'' \right) +O(h^6) $$

第三步：代入微分方程

由： $$ y''=F-ky $$ 所以： $$ y_{i+1}''=F_{i+1}-k_{i+1}y_{i+1} $$ 代入： $$ y_{i+1}-2y_i+y_{i-1} = \frac{h^2}{12} \left[ F_{i+1}+10F_i+F_{i-1} - \left( k_{i+1}y_{i+1} +10k_i y_i +k_{i-1}y_{i-1} \right) \right] $$ 把所有含 $y_{i+1}$ 的项放左边： $$ y_{i+1} + \frac{h^2}{12}k_{i+1}y_{i+1} = 2y_i-y_{i-1} + \frac{h^2}{12} (F_{i+1}+10F_i+F_{i-1}) - \frac{h^2}{12} (10k_i y_i+k_{i-1}y_{i-1}) $$ 左边： $$ \left( 1+\frac{h^2}{12}k_{i+1} \right)y_{i+1} $$ 右边整理 $y_i$ 项：

$$ 2y_i-\frac{10h^2}{12}k_i y_i = 2\left(1-\frac{5h^2}{12}k_i\right)y_i $$ 整理 $y_{i-1}$ 项： $$ -y_{i-1}-\frac{h^2}{12}k_{i-1}y_{i-1} = -\left( 1+\frac{h^2}{12}k_{i-1} \right)y_{i-1} $$ 因此： $$ \left( 1+\frac{h^2}{12}k_{i+1} \right)y_{i+1} = 2\left(1-\frac{5h^2}{12}k_i\right)y_i - \left(1+\frac{h^2}{12}k_{i-1}\right)y_{i-1} + \frac{h^2}{12} (F_{i+1}+10F_i+F_{i-1}) $$ 最终得到 Numerov 递推公式： $$ y_{i+1} = \frac{ 2\left(1-\frac{5h^2}{12}k_i\right)y_i - \left(1+\frac{h^2}{12}k_{i-1}\right)y_{i-1} + \frac{h^2}{12} (F_{i+1}+10F_i+F_{i-1}) } { 1+\frac{h^2}{12}k_{i+1} } $$ 齐次情形 $F=0$ 时： $$ y_{i+1} = \frac{ 2\left(1-\frac{5h^2}{12}k_i\right)y_i - \left(1+\frac{h^2}{12}k_{i-1}\right)y_{i-1} } { 1+\frac{h^2}{12}k_{i+1} } $$ 对二阶方程的 Taylor 展开到 6 阶，巧妙构造消去误差项：

$$y_{i+1} = \frac{2\left(1 - \frac{5h^2}{12}k_i\right)y_i - \left(1 + \frac{h^2}{12}k_{i-1}\right)y_{i-1} + \frac{h^2}{12}(F_{i+1} + 10F_i + F_{i-1})}{1 + \frac{h^2}{12}k_{i+1}}$$

💡 局部误差正比于 $h^6$，比 RK4 还高一阶！

9.3 特例：F=0 的齐次形式

当 $F(x) = 0$（齐次方程 $y'' + k(x)y = 0$）：

$$y_{i+1} = \frac{2\left(1 - \frac{5h^2}{12}k_i\right)y_i - \left(1 + \frac{h^2}{12}k_{i-1}\right)y_{i-1}}{1 + \frac{h^2}{12}k_{i+1}}$$

9.4 ⚠️ 注意事项

• 需要两个起点值 $y(0)$ 和 $y(h)$ 才能启动——不像 RK4 只要一个值 是==三点递推==
• 通常先知道 $y(0)$ 和 $y'(0)$，需用其他方法（如 Taylor 展开/RK4走一步）先求出 $y(h)$
• 比 RK4 精度高且计算量小（每步只需一次 $k(x)$ 估值）

10. Schrödinger 方程数值求解

10.1 球对称势的径向方程

$$-\frac{\hbar^2}{2m}\frac{d^2u(r)}{dr^2} + \left[V(r) + \frac{\hbar^2 l(l+1)}{2mr^2}\right]u(r) = E u(r)$$

定义有效势： $$ V_{\text{eff}}(r) = V(r)+ \frac{\hbar^2l(l+1)}{2mr^2} $$ 径向方程就是： $$ -\frac{\hbar^2}{2m}u'' + V_{\text{eff}}u = Eu $$ 整理成 Numerov 形式： $$ u''+ k(r;E)u=0 $$ 其中： $$ k(r;E) = \frac{2m}{\hbar^2} \left[ E-V_{\text{eff}}(r) \right] $$ 也就是： $$ k(r;E) = \frac{2m}{\hbar^2} \left[ E-V(r) \right] - \frac{l(l+1)}{r^2} $$ 所以 Schrödinger 径向方程天然适合 Numerov。

其中 $u(r) = rR(r)$，边界条件：$u(0) = 0,\ u(\infty) \to 0$

正好利用 上边学到的Numerov方法！！！

10.2 打靶法求解束缚态

束缚态能量要求离散！！！

单纯从 $r=0$ 往外积分有问题：

• 如果 $E$ 不准，远处通常指数发散；
• 即使 $E$ 很准，数值误差也可能激发发散解。

所以常用双向积分。

向外积分 Outward

从小 $r=\epsilon$ 开始，使用原点行为： $$ u_{\text{out}}(r)\sim r^{l+1} $$ 可以取： $$ u(\epsilon)=\epsilon^{l+1} $$ 然后用 Taylor 或 RK4/Numerov 得到第二个点，开始向外积分。

向内积分 Inward

在足够大的 $R_{\max}$ 处，束缚态远处衰减。

如果远处势能趋于 0，且束缚态 $E<0$，则： $$ u''+\frac{2mE}{\hbar^2}u\approx 0 $$ 因为 $E<0$，令： $$ \kappa=\sqrt{\frac{-2mE}{\hbar^2}} $$ 则： $$ u''-\kappa^2 u=0 $$ 解为： $$ u(r)=Ae^{-\kappa r}+Be^{\kappa r} $$ 束缚态要求不发散，所以： $$ u(r)\sim e^{-\kappa r} $$ 因此在 $R_{\max}$ 取： $$ u_{\text{in}}(R_{\max})=e^{-\kappa R_{\max}} $$ 然后从 $R_{\max}$ 向内积分到匹配点 $r_m$。

匹配条件

因为线性方程的解可以整体乘任意常数，所以不能直接要求： $$ u_{\text{out}}(r_m)=u_{\text{in}}(r_m) $$ 因为振幅可以缩放。

真正关键是形状要一致，也就是对数导数一致： $$ \frac{u_{\text{out}}'(r_m)}{u_{\text{out}}(r_m)} = \frac{u_{\text{in}}'(r_m)}{u_{\text{in}}(r_m)} $$ 定义匹配函数： $$ F(E)= \frac{u_{\text{out}}'(r_m)}{u_{\text{out}}(r_m)} - \frac{u_{\text{in}}'(r_m)}{u_{\text{in}}(r_m)} $$ 本征能量满足： $$ F(E)=0 $$ 等价地，可以写成 Wronskian： $$ W(E) = u_{\text{in}}u_{\text{out}}' - u_{\text{out}}u_{\text{in}}' $$ 如果两者对数导数相等，则： $$ u_{\text{in}}u_{\text{out}}' = u_{\text{out}}u_{\text{in}}' $$ 所以： $$ W(E)=0 $$ 因此求本征能量就是求： $$ W(E)=0 $$ 或者： $$ F(E)=0 $$ 可以用：

• 二分法；
• 割线法；
• 牛顿法；
• 节点数扫描。

寻找束缚态本征解流程：
1. 从 r=0 向外积分 (outward integration)
2. 从 r=R∞ 向内积分 (inward integration)  
3. 在大距离 ρₘ 处匹配 (matching point)
4. 要求两个解的 Wronskian 满足：u_in·u_out' - u_out·u_in' 连续
5. 调节能量 E，重复上述步骤
→ 找到使匹配条件成立的 E = 本征能量

flowchart LR
    E[猜测能量 E] --> Out[0→匹配点<br/>向外积分]
    E --> In[Rmax→匹配点<br/>向内积分]
    Out --> Match[匹配点处比较<br/>对数导数]
    In --> Match
    Match --> Check{连续?}
    Check -->|否| Adjust[调节 E] --> E
    Check -->|是| Done[✓ E 为本征值]

💡 物理直觉：束缚态能量是离散的——就像吉他弦只有特定频率才能共振。向外积分和向内积分只有在正确的能量 $E_n$ 下才能光滑连接。

在一维或径向 Sturm-Liouville 问题中，束缚态有节点定理：

• 基态没有节点；
• 第一激发态有一个节点；
• 第二激发态有两个节点；
• 依此类推。

所以扫描能量时，可以通过节点数判断你找到的是哪一个本征态。

10.3 散射态与共振态

• 散射边界条件：$u(r) \sim \sin(kr + \delta)$，提取相移 $\delta$
• 共振态：能量为复数 $E - i\Gamma/2$，相移突然增加 $\pi/2$
• 求解方法：复标度法 (Complex-scaling)、稳定化法 (Stabilization)、S-矩阵极点法

11. Green 函数与 Wronskian

11.1 核心问题

Green 函数是线性微分方程的“响应函数”。

考虑： $$ L[y]=f(x) $$ 如果能找到 $G(x,\xi)$ 满足： $$ L_x[G(x,\xi)]=\delta(x-\xi) $$ 并满足同样的边界条件，那么解就是： $$ y(x)=\int_a^b G(x,\xi)f(\xi)\,d\xi $$ 设： $$ L[y]=y''+k(x)y $$ 边界条件： $$ y(a)=0,\qquad y(b)=0 $$ 我们要构造： $$ G''(x,\xi)+k(x)G(x,\xi)=\delta(x-\xi) $$ 当 $x\neq \xi$ 时： $$ \delta(x-\xi)=0 $$ 所以 Green 函数在左右两边都满足齐次方程： $$ G''+kG=0 $$

11.2 构造左右基本解

令 $u_<(x)$ 是满足左边界条件的齐次解： $$ u_<(a)=0 $$ 令 $u_>(x)$ 是满足右边界条件的齐次解： $$ u_>(b)=0 $$ 于是 Green 函数可以写成： $$ G(x,\xi) = \begin{cases} A(\xi)u_<(x), & x<\xi\\ B(\xi)u_>(x), & x>\xi \end{cases} $$ 因为：

• $x<\xi$ 时要满足左端边界；
• $x>\xi$ 时要满足右端边界。

11.3 连续条件

Green 函数本身必须连续。

如果 $G$ 在 $x=\xi$ 处跳跃，那么 $G'$ 会有 delta，$G''$ 会有 delta 的导数 $\delta'$，而方程右边只有 $\delta$，没有 $\delta'$。

所以： $$ G(\xi^-,\xi)=G(\xi^+,\xi) $$ 即： $$ A u_<(\xi)=B u_>(\xi) $$

11.4 导数跳跃条件

对方程在 $\xi-\epsilon$ 到 $\xi+\epsilon$ 积分： $$ \int_{\xi-\epsilon}^{\xi+\epsilon} \left[ G''+kG \right]dx = \int_{\xi-\epsilon}^{\xi+\epsilon} \delta(x-\xi)\,dx $$ 右边： $$ 1 $$ 左边第一项： $$ \int G''dx = G'(\xi^+)-G'(\xi^-) $$ 第二项： $$ \int kG\,dx \to 0 $$ 因为区间长度趋于 0，且 $G$ 有限。

所以： $$ G'(\xi^+)-G'(\xi^-)=1 $$ 代入左右表达式： $$ G'(\xi^-)=A u_<'(\xi) $$ 因此： $$ B u_>'(\xi)-A u_<'(\xi)=1 $$

11.5 解出 $A,B$

我们有： $$ A u_<(\xi)=B u_>(\xi) $$ 所以： $$ A=\frac{B u_>(\xi)}{u_<(\xi)} $$ 代入跳跃条件： $$ B u_>'(\xi) - \frac{B u_>(\xi)}{u_<(\xi)}u_<'(\xi) = 1 $$ 提取 $B$： $$ B \left[ u_>'(\xi) - \frac{u_>(\xi)u_<'(\xi)}{u_<(\xi)} \right] = 1 $$ 通分： $$ B \frac{ u_<(\xi)u_>'(\xi) - u_>(\xi)u_<'(\xi) } {u_<(\xi)} = 1 $$ 定义 Wronskian： $$ W(\xi) = u_<(\xi)u_>'(\xi) - u_>(\xi)u_<'(\xi) $$ 于是： $$ B\frac{W(\xi)}{u_<(\xi)}=1 $$ 所以： $$ B=\frac{u_<(\xi)}{W(\xi)} $$ 再由连续条件： $$ A=\frac{u_>(\xi)}{W(\xi)} $$ 因此： $$ G(x,\xi) = \begin{cases} \dfrac{u_<(x)u_>(\xi)}{W(\xi)}, & x<\xi\\[8pt] \dfrac{u_<(\xi)u_>(x)}{W(\xi)}, & x>\xi \end{cases} $$ 如果方程没有一阶导数项，Wronskian 是常数，所以也可以写成： $$ W $$ 而不写 $W(\xi)$。

11.7 Green 函数和本征值的关系

如果： $$ W=0 $$ 说明 u_{<} 和 u_{>} 线性相关。

也就是说存在非零解同时满足左右边界条件。

这正是本征值问题的条件。

所以：

Green 函数的分母 $W$ 在本征值处为零。

这也解释了物理中常说：

Green 函数在能量本征值处有极点。

构造流程（配合 Numerov）：

1. 从 a→b 向外积分齐次方程，得到 u<(x)|满足 a 端边界条件
2. 从 b→a 向内积分齐次方程，得到 u>(x)|满足 b 端边界条件
3. 计算 Wronskian W = u<·u>' - u>·u<'
4. 构造 G(x,ξ)，积分求最终解

💡 物理直觉：Green 函数 = 线性响应函数。$G(x,\xi)$ 是"$\xi$ 处的单位脉冲在 $x$ 处引起多大的响应"，与外场 $f$ 无关。这是整个线性响应理论的数学基础。

12. 有限差分法（Finite Difference）

12.1 基本思想

将导数用差分代替，把微分方程变成代数方程组。

一阶和二阶导数的中心差分近似（精度 $O(h^2)$）：

$u'(x) = \frac{u(x+h)-u(x-h)}{2h} +O(h^2)$

$u_i'' \approx \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$

12.2 📐 数值例子：一维泊松方程

$$-u''(x) = f(x), \quad u(0)=a,\ u(L)=b$$

离散化【$N$ 个内点（未知），$h = L/(N+1)$ 】：

通常写成： $$ \frac{-u_{i-1}+2u_i-u_{i+1}}{h^2}=f_i $$ 其中： $$ f_i=f(x_i) $$

内部点方程

对于 $i=1$： $$ \frac{-u_0+2u_1-u_2}{h^2}=f_1 $$ 因为： $$ u_0=a $$ 所以： $$ \frac{2u_1-u_2}{h^2} = f_1+\frac{a}{h^2} $$ 对于中间点 $i=2,\dots,N-1$： $$ \frac{-u_{i-1}+2u_i-u_{i+1}}{h^2}=f_i $$ 对于 $i=N$： $$ \frac{-u_{N-1}+2u_N-u_{N+1}}{h^2}=f_N $$ 因为： $$ u_{N+1}=b $$ 所以： $$ \frac{-u_{N-1}+2u_N}{h^2} = f_N+\frac{b}{h^2} $$

对于内部点转化为线性方程组 $A\mathbf{u} = \mathbf{b}$： $$ \frac{1}{h^2} \begin{pmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \cdots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 \\ 0 & 0 & 0 & -1 & 2 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_N \end{pmatrix} = \begin{pmatrix} f_1 + a/h^2 \\ f_2 \\ f_3 \\ \vdots \\ f_N + b/h^2 \end{pmatrix} $$

说人话：一个二阶 ODE 的边值问题变成了一个三对角线性方程组！用 Thomas 算法 $O(N)$ 秒解，干净利落。

import numpy as np

def solve_poisson_1d(a, b, f, L, N):
    """解 -u'' = f, u(0)=a, u(L)=b"""
    h = L / (N + 1)
    x = np.linspace(0, L, N+2)
    
    # 构造三对角矩阵 + RHS
    diag = 2 * np.ones(N) / h**2
    off_diag = -np.ones(N-1) / h**2
    rhs = f(x[1:-1])
    rhs[0] += a / h**2
    rhs[-1] += b / h**2
    
    # Thomas 算法 (三对角追赶法)
    # 向前消元
    c_prime = np.zeros(N-1)
    d_prime = np.zeros(N)
    c_prime[0] = off_diag[0] / diag[0]
    d_prime[0] = rhs[0] / diag[0]
    for i in range(1, N):
        if i < N-1:
            denom = diag[i] - off_diag[i-1] * c_prime[i-1]
            c_prime[i] = off_diag[i] / denom
        else:
            denom = diag[i] - off_diag[i-1] * c_prime[i-1]
        d_prime[i] = (rhs[i] - off_diag[i-1] * d_prime[i-1]) / denom
    
    # 回代
    u_interior = np.zeros(N)
    u_interior[-1] = d_prime[-1]
    for i in range(N-2, -1, -1):
        u_interior[i] = d_prime[i] - c_prime[i] * u_interior[i+1]
    
    u = np.zeros(N+2)
    u[0], u[-1] = a, b
    u[1:-1] = u_interior
    return x, u

# 示例：f(x)=sin(x), a=0, b=0, L=π
# u_exact(x) = sin(x) → 验证数值解

12.3 有限差分法 vs 打靶法

💡 口诀：

一维光滑问题，打靶法很灵活； 线性边值问题，有限差分很稳； 多维边值/PDE，有限差分、有限元更主流。

13. 应用实例：TOV 方程与中子星

13.1 TOV 方程 (Tolman-Oppenheimer-Volkoff)

广义相对论中描述中子星结构的耦合 ODE 组：

$$\begin{aligned} \frac{dP}{dr} &= -\frac{G}{r^2}\frac{[\rho(r) + P(r)/c^2][m(r) + 4\pi r^3 P(r)/c^2]}{1 - 2Gm(r)/(rc^2)} \\ \frac{dm}{dr} &= 4\pi r^2 \rho(r) \end{aligned}$$

13.2 状态方程 (EOS)

核物理提供密度-压强关系。简化多方形式： $$P = K \rho^\gamma$$

典型参数：$\gamma = 2.75$，$K = 1.982 \times 10^{-6}$ [cgs]

13.3 求解策略

初始条件 (r=0):
  ρ(0) 给定（中心密度）
  P(0) = P(ρ(0))  从 EOS 得到
  m(0) = 0

用 RK4 从 r=0 向外积分
→ 每次步进：r → r + dr
  更新 P(r), m(r)（耦合方程组）
→ 直到 P(R) = 0（压力降到零）
→ 此时 R = 中子星半径，M = m(R) = 中子星质量

改变中心密度 ρ(0)，重复 → 得到质量-半径关系
→ 最大质量 ≈ 2-3 倍太阳质量（取决于 EOS）

💡 物理直觉：用 RK4 从恒星中心往外积，密度（因此压力）一路下降。到表面压力为 0 → 这就是星体边界。不同的中心密度 → 不同的总质量 → 找到最大质量（超过此值星体会坍缩为黑洞）。

13.4 前沿研究

• Chong-Ji Jiang et al., Chin. Phys. Lett. 2021, 38(5): 052101
• 中子星物质状态方程仍是核物理前沿
• 引力波观测 (GW170817) 给 EOS 带来新的约束

14. 总结对比表

14.1 初值问题方法总览

14.2 边值问题方法总览

14.3 关键联系与口诀

flowchart TD
    subgraph 初值问题
        E[Euler] --> T[Taylor 展开]
        T --> RK[Runge-Kutta]
        RK --> PC[Predictor-Corrector]
        PC --> TDDFT[TDDFT 等长时间计算]
    end
    subgraph 边值问题
        ST[打靶法 Shooting] --> IVP2[转化为初值问题]
        IVP2 --> RK2[用 RK4 积分]
        Numerov2[Numerov 方法] --> Schro[薛定谔方程]
        Green[Green 函数] --> Numerov3[Numerov 积分]
        FD[有限差分法] --> Matrix[线性方程组]
    end
    subgraph 特殊问题
        Stiff[Stiff 方程] --> Implicit[隐式方法<br/>梯形法/Gear]
    end

📚 核心记忆卡片

Akuiro 整理补充 · 2026-07-01