板书小抄

“rtol = 1e−6 的这条和参考解，哪一条是对的？”

停一下等回答。脱钩之前参考解对；脱钩之后两条都只是吸引子上的某一条轨道，“哪条对”这个问题本身没了意义。后两个间隔都是 8 左右：每收紧 1000 倍只多 8，这个 8 第 ② 部分算。目录页只请大家记住 2.55。

容差收紧十亿倍，可信时间只从 5 走到 27。

“这组方程从哪来，怎么放进计算机？”

锅里的水底下烧：热的变轻想上去，冷的想下来；粘性 ν 拖住流动，热扩散 κ 抹平温差。Ra 小拦得住，水不动；过了临界值长出一排对流卷。能拧的只有 ΔT，Ra ∝ ΔT，这个旋钮后面叫 ρ。

θ 是实际温度减掉传导直线：水不动 θ = 0，所以 θ 量的是对流本身。

两条 PDE，未知量是函数，要先变成有限个数。

线方法：先空间离散、再时间推进，和声学 FDTD 一样；只是第 1 步用 Fourier 模态，只留 3 个，N = 3，小到每一步都看得清。

大写 X、Z 是坐标，小写 x、y、z 是只随时间变的振幅（H = 1，k = πa）。z 的模态只随竖直 Z 变：中段被搅匀拉平、近板变陡，所以 z 就是对流多运走的热。

x 是转速，y 是温差，z 是多运的热。

“只留三项，为什么能自己闭合？扔掉了什么？”

② 是 z 非留不可的理由：x、y 一乘直接生出 z。不留 z，方程只剩线性部分，ρ > 1 时振幅指数增长，没有东西让它停；z 是非线性的刹车。

③ 是截断发生的唯一地方。被问到可现场验 ②：积化和差后含 cos 2kX 的项正好抵消，只剩 sin 2πZ。

截断就是扔掉 sin 3πZ；z 是非线性的刹车。

刚过 ρ = 1，xz 是小量乘小量，扔掉的那块也小，截断有渐近依据；到 28，扔掉的强迫和驱动 ρx 同一量级，真实流体里 sin 3πZ 一定长起来。传热是第二个迹象：截断把 Nu 封顶在 3。

第 10 页一句带过：线性项只有驱动和扩散，非线性两项 −xz、xy 都是 u·∇θ，混沌全来自搬运。

恰好就是这组方程的装置：Malkus 水车、热虹吸环、单模激光（Haken 1975）。

后面的结论关于方程本身，不关于天气。

“为什么是 8/3 这么怪的数？”

8/3 是最先失稳的那种卷的几何；σ 和水同量级（室温约 7），还得大于 β + 1，否则 C± 永不失稳；28 刚过 24.74。量级（被问再说）：1 cm 水层上下差 1 K，Ra 约 1.5×10⁴；一个时间单位约一分钟；刚性板临界值是 1708。

8/3 是几何，10 是水，28 刚过失稳。

“怎么让计算机往前走一步？”

左边代码就是全部计算：一个 f、一个循环。全篇实算只用 DOP853，只拧 rtol，开场演示就在拧它。过渡：方程有了，积分器也有了；谈误差之前，先看方程自己会做什么。

每一步都很准；问题是误差怎么累积。

每拖之前先猜：“ρ 小于 1 会怎样？刚过 1？一直加热下去？”

0.5 不动；过 1 卷自己长出来，转向随机；10 定常卷更强；过 24.74，x 在正负间跳，停不下来。24.5 推一下（z 加 15）就不再回来：两种结局并存，看初值。

拧 ρ，从静止到定常卷，再到停不下来的翻转。

“刚才停下来的点，能不能不积分，直接从方程里找出来？”

O 是静止导热，C± 是定常卷，正负是转向；(x, y) → (−x, −y) 对称，所以成对。原点附近 xz、xy 是小量乘小量，扔掉就是线性化，z 自己分出去。ρ = Ra/Ra_c 这样归一化，就是为了让阈值落在 1。

ρ = 1，原点从稳定变成鞍点，C± 诞生。

对着第 16 页四格读：0.5 一个点；10 原点空心、两条轨道螺旋落进 C±；24.5 一条在两叶间游荡、一条困在 C+ 旁；28 三个点全空心。

第 17 页分岔图的纵轴：固定 ρ，把 z 每个波峰竖着点上去。ρ_H 以下那一列收到 ρ − 1，正是 C± 的高度：虚部让轨道绕着转出峰，负实部让峰逐圈收拢。带在 ρ ≈ 23.8 就出现（23.8 到 24.06 是长暂态），比 ρ_H 早。

“一直加热，定常卷什么时候站不住？”

系数全正挡住的只是实根：反例 λ³ + λ² + λ + 2，实根约 −1.35，三根之和 −1，复根实部 +0.18。根随 ρ 连续地动，到右半平面必须过虚轴；a₀ > 0，λ = 0 不是根，只能在 ±iω 成对穿。分母 σ − β − 1 说明为什么要 σ > β + 1。

亚临界：穿越时没有小极限环接住，轨道一离开 C± 就跳到远处的吸引子。

第一条：停不下，ρ = 28 时三个不动点全不稳。

三个系数就够验算，不用解三次方程。Vieta：a₂ = −(根之和) = −tr J，a₀ = −(根之积) = −det J。

ρ 每加 1，a₀ 涨 160/3 ≈ 53，a₂a₁ 只涨 41/3 × 8/3 ≈ 36，所以一过 ρ_H 就一直是 a₀ 大。

两个不动点的特征值加起来都是 −13.67。先记下这个数：第 23 页它是 tr J。

“三个不动点都站不住，轨道在两叶之间怎么换？”

竖一块平面，轨道每从同一侧穿过一次记一个点，其余全扔：三维的流变成平面上的映射，时间换成“第几圈”。每一叶是一张薄片，切过平面是一条细线，细到看不出宽度，就是 D_KY = 2.06。

一个点是一步，不是一个位置。0→1→2 沿绿色往外爬到门槛；3 是 13.3，只越过 0.06，被甩到对面最外圈 −15.3；4 越过 2，反而落到右叶中心附近的 11.4。截面降到一维就不受“流线不交叉”的限制，点一下一下跳，映射上可以有混沌。

近邻每圈乘 |F′| > 1；“乘多少”变成一个数就是 Lyapunov 指数，要算它先在任意一点线性化（下一页）。

三维的流压成一条一维映射，映射处处拉伸。

“轨道上任意一点旁边差一点点，这一点点下一刻怎么变？”

J 的第 i 行回答：第 i 个变量的变化率对 x、y、z 各多敏感。线性项求导得常数，xz、xy 求导留下 x、y、z，所以轨道走到哪 J 就是哪个矩阵。

J 里有 x、y、z，所以它随轨道上的点变。

“轨道走到这一点，扰动在这里长得多快？”

|δu|² 求导得 2δu·Jδu，除以 2|δu|² 就是 e·Je。只给一个点不够，还要说扰动指向哪。特征值说“一直乘同一个矩阵，最后谁赢”；e·Je 说“此刻这个方向怎样”。下一刻矩阵换了，本征方向也换了。

窗口 τ 上平均 = 有限时间指数：τ = 0.5 标准差 2.60，43% 的窗口在缩；τ → ∞ 收拢到 λ₁。图上红色在两叶交界与外缘。

局部指数随点、随方向跳，只有长期平均才是 λ₁。

“J 在不动点上看特征值就行，沿着轨道为什么不行？”

每走 Δt 乘一次 I + J(u)Δt，沿轨道是一长串各不相同的矩阵，没有共用的本征方向。随手跟三个向量，全倒向 λ₁：第 22 页不做 QR，t = 30 时夹角 3×10⁻¹⁴（心算 e^{−(λ₁−λ₂)t} = e⁻³⁰ ≈ 10⁻¹³），三条斜率 0.98 / 1.01 / 1.05 全是 λ₁；做 QR 后分开成 0.92 / −0.003 / −14.58。

为什么每 Δt = 1 就做：夹角到 10⁻¹⁶ 就丢，ln(10¹⁶)/0.9 ≈ 41，和双精度预报极限是同一个数；长度到 e⁷⁰⁹ 溢出。

第二条：一个方向在拉伸，λ₁ ≈ 0.9。

“三个指数加起来是多少？能不能不算？”

换一组正交坐标轴不改迹。Lorenz 特殊在对角线 −σ、−1、−β 是粘性和两项热扩散，没有变量：搬运只挪体积，耗散才缩体积。和式守的是变分方程、积分精度、QR 记账；三个数怎么分它管不到，λ₁ 单独的依据是第 26 页系综斜率 0.875。

回扣那一行是 Abramov 公式：截面映射的指数 = λ₁ × 平均回归时间；拿 fig14 的截面数据拟合局部斜率核过，平均 ln|F′| = 0.68。停不下、跑不掉（备用 2）、体积为零 → 奇怪吸引子。

第三条：体积缩到零，Σλ = tr J = −13.67。

“算得更准，能多算多久？”

第 25 页先跑 ε = 1e−6；再问“缩小一百万倍到 1e−12，能多看多久？”猜完再扫：每条线平移，6 × 2.55 ≈ 15，不是一百万倍的时间，是加 15。

第 26 页：16 个吸引子上的状态各推 4 个方向 1e−9，共 64 对；读三段：对齐、直线、在 20 左右走平。斜率 0.875 对 QR 的 0.904，差 3%，没有共用代码。第 27 页：rtol 也是 δ；拟合斜率 1.131 对 1/λ₁ = 1.107，差 2%（有效初始误差与 rtol 同阶，不严格成正比）。

老实交代：单个起点相邻间隔在 0.2 到 7 之间跳，所以要多起点取中位。

精度乘 10，时间只加 2.55。

rtol 和双精度是两回事。就算算法完美，初值只能存到 1e−16，t > 41 算的就不是那个初值的轨道：“哪条对”在双精度下没有答案。到这里全是坏消息。

t > 41，单条轨道在双精度下没有答案。

第 28 页算账：8 阶方法步数随 rtol 的 −1/8 次方涨，8 个数量级正好 10 倍。10 倍的计算换来 5.7 到 26.9；再要 20 个时间单位得再收紧 8 个数量级，双精度已经到头。

第 38 页 Kepler（e = 0.6，RK4 步长取 Verlet 的 4 倍，同等代价）：Verlet 能量误差 7.4×10⁻⁴ 跑 1000 周期不涨，RK4 线性漂移，第 33 周期超过；可信时间对代价是幂律，斜率 2.0 / 2.5。

精度换时间：守恒是乘法，混沌是加法。

嵌入式：同一组 k_i 两套权重，两个解之差就是本步误差估计，几乎不花额外 f；err ≤ 1 接受，下一步 h ← 0.9 h err^−1/p。

被问“为什么用 DOP853”：rtol = 1e−4 时两者 f 次数持平，容差越紧 DOP853 越划算。刚性时步长被稳定域卡住，s / 交点三者差不多，阶数换不来步数。

积分器决定代价，λ₁ 决定 2.55。

增长模（备用 6 的红点）永远在稳定域外，那是格式在正确复现真实的增长；要求全谱 |R| ≤ 1 只会得到 h = 0。稳定域只约束衰减模。

衰减模的振幅早衰减没了，误差估计看不见它，控制器不管它的精度，把它带出精度域，但还在稳定域里：步长由慢模的精度定。

刚性反例（备用 9）：y′ = −1000(y − cos t)，T = 10，步数 ≈ 10 × 1000 / 交点：RK45 3020（实测 3025），DOP853 1560（实测 1567），与 rtol、阶数都无关；隐式 Radau 只要 112 次。

步长由慢模的精度定，衰减模只要还在稳定域里就没事。

被问“步长再小一点不就行了？”（备用 4）：数值稳定性是方法加步长的性质，步长能改；Lyapunov 指数是流的性质，步长改不了。步长变小只减少注入的误差量，这些误差随后照样按 e^λt 放大。

被问“红色区域减步长有没有用？”：有用，理由是少往那里注入误差，不是压住 λ。

“单条轨道不可信，那还有什么可信？”

第 31 页多停：左边问“t 时刻你在哪”，取决于初值，被 e^λ₁t 放大；右边问“从这么高的峰出发，下一个峰多高”，是方程本身的性质。映射斜率绝对值处处大于 1：左支 +1.1 到 +2.6，右支 −1.15 到 −4.4。

同一组数据，按时间看是废品，按映射看完好。

“换一类系统，答案还一样吗？”

第 36 页心算：ρ = 350 周期 0.39，λ₂ = −0.73，每周期乘 e^−0.28 ≈ 0.75；ω = 2/5 每次转 144°，5 步闭合；ω = 0.618 每次转 222.5°（反着 137.5°，黄金角），永不闭合，就是李萨如图。

第 37 页：守恒量沿轨道是常数且连续；吸引盆里每条轨道都贴到吸引子上，混沌轨道走遍吸引子，于是它在整个吸引盆上是常数，什么也区分不了。

守恒量记得初值，吸引子让轨道忘掉初值。

第 38 页 (c) 老实讲：同等代价下 RK4 的线在 Verlet 上面；辛格式买到的不是更准，而是不漂，跑得越久越占优。

第 40 页把黑板擦一遍：一，没有守恒量，体积以 −13.67 缩；二，精度乘 10 只多 2.55，天花板 41；三，Lorenz 映射、不变密度不随 rtol 变。

守恒系统看不变量，混沌系统看统计量；精度换时间，前者乘法，后者加法。