从对流方程
到可信的数值解

底下的流体受热变轻想上浮，上面的冷流体想下沉。粘性和热扩散会抵抗这种翻动。

两股力量之比叫 Rayleigh 数：

• Ra 小于临界值 Ra_c：流体静止，热量只靠传导
• Ra 大于 Ra_c：静止失稳，长出一排对流卷
• 另一个无量纲数 Pr = ν/κ：动量与热量各自扩散的快慢之比

Saltzman（1962）、Lorenz（1963）只看一个方向上的二维卷。

Boussinesq 近似：密度只在浮力项里随温度变，其余地方当常数。

• 只看 X–Z 平面里的二维卷：令 \(u = -\partial_Z\psi,\ w = \partial_X\psi\)，质量守恒自动满足
• 对动量方程取旋度，压力 p 消掉
• 温度拆成两份：\(T = \underbrace{T_0 + \Delta T\,(1 - Z/H)}_{\text{静止时的传导直线}} + \theta\)

剩下两个未知场：流函数 ψ(X, Z, t) 管流动，θ(X, Z, t) 是温度偏离传导直线的那部分。流体静止时 θ ≡ 0。

\(J(a,b) = \partial_X a\,\partial_Z b - \partial_Z a\,\partial_X b\)，就是 \(\mathbf u\cdot\nabla b\)。这是 Lorenz 1963 的式 (17)–(18)。未知量是两个函数，所以是无穷维的。

后面所有结论都关于这个数学对象，不涉及流体或天气的主张。

但有几个物理系统，化简后恰好就是这组方程：

• Malkus 水车：x 是转速，y、z 是桶中水量分布的一阶 Fourier 系数
• 热虹吸环：竖直圆环管里的对流，x 是环流速度
• 单模激光（Haken 1975）：x 是光场，y 是极化，z 是粒子数反转

def f(u):                # 就是那三行方程
    x, y, z = u
    return np.array([σ*(y-x), x*(ρ-z) - y, x*y - β*z])

u = u0                   # 显式 Euler：一阶
for n in range(N):
    u = u + h * f(u)

def rk4_step(u, h):      # 经典 RK4：四阶
    k1 = f(u)
    k2 = f(u + h/2 * k1)
    k3 = f(u + h/2 * k2)
    k4 = f(u + h * k3)
    return u + h/6 * (k1 + 2*k2 + 2*k3 + k4)

仓库里的 integrate.rk4_step 就是这几行（多一个时间参数）。后面的计算和现场演示都用 DOP853，只调 rtol。

线性项求导得常数；xz、xy 求导留下 x、y、z，J 随轨道上的点变。

• 不动点上　J 是常数，看特征值：不动点与 Hopf 两页已经用过
• 沿轨道　J 每一点都不同，没有固定的本征方向：局部指数、边算边正交化
• 对角线　tr J = −(σ+1+β)，与点无关：三个指数之和

直觉：局部误差很小的数值轨道，虽然不是你要的那条真解，可能紧跟着另一条初值稍有不同的真解。

那么画出来的蝴蝶是真的蝴蝶，只是不是"你那只"。

对一致双曲系统，这是定理（影随引理）。

影随只提示方向。真正的辩护必须可测量。

表的第一行

无耗散波动方程是哈密顿系统。Yee 的 FDTD 格式里 E 与 H 在时间上交错半格，那就是蛙跳格式，也就是辛积分器。

FDTD 能长时间运行而不出现数值能量爆炸或人工衰减，原因在结构。

表的第三行

非线性声学、声致湍流一旦进入混沌，单次实现的波形不该当作可复现的预报。

该报告的是谱、相关函数、统计矩。

• Lorenz 没有守恒量，体积以 tr J = −13.667 坍缩。
• 可预报时间被 λ₁ = 0.904 锁死在约 41 个时间单位，精度只换来对数的时间。
• 不变密度、Lorenz 映射对 rtol 不敏感。

全部数字与图都可复现，关键断言写成了测试：

uv sync --all-groups
uv run pytest
uv run lorenz63 info
uv run lorenz63 figure all

RK45 即 scipy 的 Dormand–Prince 5(4)；DOP853 是 Hairer–Wanner 的 8(5,3)。演示里的 DOP853 由 scipy 逐行移植。

步数 \(\propto \text{rtol}^{-1/p}\)。rtol 收紧 10⁶ 倍：

积分器决定代价；每 10 倍 rtol 换 2.55 个时间单位，由 λ₁ 决定。

方括号在椭球 \(\rho x^2 + y^2 + \beta(z-\rho)^2 = \beta\rho^2\) 之外为正，V 在那里单调下降。

所有轨道最终被困在一个有界区域里。

线性化 \(\dot{\delta u} = J\,\delta u\)，对角化后每个方向都是 \(\dot y = \lambda y\)。

\(R_\text{RK4}(z) = 1 + z + \tfrac{z^2}{2} + \tfrac{z^3}{6} + \tfrac{z^4}{24}\)，是 \(e^z\) 的截断。

h·Re λ：一步里幅度的对数变化。h·Im λ：一步转过的弧度（C±：0.10，62 步一圈）。

最快的衰减模 λ = −21.85，时间尺度 0.046：

R(z) 由 Butcher 表算出；上限 = 交点 / 21.85。与 rtol 无关。

DOP853 在吸引子上的最大步长，占稳定上限：

全文用的 10⁻⁴ 到 10⁻¹³ 都在稳定域里面：步长按精度选。

常见做法：画"误差 vs 步长"，期待一条 O(h^p) 直线。

超过可预报时间后误差饱和到吸引子直径，图会变成一条水平线，看起来像格式失去了精度阶。

test_integrate.py::
test_tighter_tolerance_is_more_accurate_before_saturation