3.1 3.1 ODE 建模与数值解：反应釜降温工单

3.1 ODE 建模与数值解：反应釜降温工单

本节摘要：常微分方程建模的套路是"找出状态量、写出每个状态的变化率由什么决定"。本节以间歇反应釜的夹套降温为工单，列状态方程、手写欧拉法看清数值方法的骨架、再切换到龙格-库塔与自适应步长求解器，最后讨论刚性方程为什么需要专门的求解器。

学习目标

阅读完本节，你应当能够：

1. 把连续变化的物理过程翻译成状态方程组；
2. 手工实现欧拉法与改进欧拉法并解释误差阶;
3. 判断问题是刚性的还是常规的，选对求解器。

工单：给反应釜降温

一批放热反应结束后，工艺要求两小时内把釜内温度从 95 摄氏度降到 40 度以下，夹套通冷却水。温度怎么随时间变化？能不能赶上交接班？这里的状态量只有一个：釜内温度 T。它的变化率由能量平衡决定：釜内热量流失速率正比于釜温与冷却水温之差。

写成方程：温度对时间的导数等于负的传热系数参数乘以（T 减去冷却水温）。这个一阶线性方程其实有解析解——指数衰减。先别高兴：一旦冷却水温随时间变化、或反应余热还在释放，解析解立刻消失，数值方法是唯一出路。而实际工单几乎总有这两项。

数值解法的基本想法朴素到不像数学：知道当前时刻的温度和变化率，就用"当前位置加斜率乘步长"走到下一时刻。这就是欧拉法。它用折线逼近曲线，误差与步长成正比——步长减半误差减半，看起来能用，但要四位有效数字就得把步长压到很小。

import numpy as np

def cooling(t, T, k=0.03, Tc=20.0):
    """釜温变化率：正比于与冷却水的温差"""
    return -k * (T - Tc)

def euler(f, T0, t0, t1, h):
    """显式欧拉法：折线逼近"""
    ts = np.arange(t0, t1 + h/2, h)
    Ts = np.empty_like(ts)
    Ts[0] = T0
    for i in range(len(ts) - 1):
        Ts[i+1] = Ts[i] + h * f(ts[i], Ts[i])
    return ts, Ts

T_exact = lambda t: 20 + 75 * np.exp(-0.03 * t)   # 解析解
t_end = 120

for h in [8, 4, 2, 1]:
    ts, Ts = euler(cooling, 95.0, 0, t_end, h)
    err = abs(Ts[-1] - T_exact(t_end))
    print(f"步长 {h:>3}: 终温数值 {Ts[-1]:.4f}, 误差 {err:.4f}")
​

运行后能看到误差大致随步长线性缩小——一阶方法的标志。改进办法是让折线的斜率取"起点与终点斜率的平均"（改进欧拉，二阶），或更精细的四点加权（经典龙格-库塔，四阶）。四阶方法步长减半误差降到十六分之一，这是它统治教科书八十年的原因。

自适应步长与库的正确用法

真实求解不自己造轮子。科学计算库的求解器带自适应步长：误差小时放大步长赶进度，误差大时缩小步长保精度。用法上有一个新手高频翻车点——参数顺序，以及刚性判断。

import numpy as np
from scipy.integrate import solve_ivp

def cooling(t, T, k=0.03, Tc=20.0):
    return -k * (T - Tc)

sol = solve_ivp(cooling, [0, 120], [95.0],
                method='RK45', rtol=1e-8, atol=1e-10,
                dense_output=True)

t_check = np.array([30, 60, 120])
T_num = sol.sol(t_check)[0]
T_ref = 20 + 75 * np.exp(-0.03 * t_check)
print("数值解:", np.round(T_num, 6))
print("解析解:", np.round(T_ref, 6))
print("内部步数:", len(sol.t))
​

容差设到 10 的负 8 次方后，数值解与解析解在小数点后六位一致，而求解器内部只用了几十步——自适应的优势是"该密则密、该疏则疏"。

刚性是另一类病情：系统里同时存在快变量与慢变量。典型如带强放热余热的反应釜——釜温降得慢，夹套壁温响应快，两者的时间尺度差两个数量级。此时常规显式方法被迫用极小步长去稳定那个快变量，慢变量明明可以大步走。刚性问题要用隐式方法（如库里的 Radau 或 BDF）：每步解一个小型方程组，换来的是大步长下的稳定性。判断办法很土：先跑 RK45，若内部步数远超预期、步长被压得极小，换 BDF 再试，耗时骤降即确认刚性。

⚠️ 常见坑：把求解区间拆成循环手动调用求解器来"控制步长"。自适应求解器本身就是干这个的，手动拆区间既破坏误差估计又拖慢速度。要输出等间隔结果，用稠密输出（dense_output）事后插值。

💡 关键直觉：欧拉法的几何图像值得刻在脑子里——所有高阶方法都是"如何更聪明地选斜率"的变体。龙格-库塔在区间内多采样几个斜率做加权平均，误差阶就上去了。理解这一点，库的参数（容差、方法名）就不再是黑盒开关。

求解路线决策

建模检查单与常见疏漏

收尾给一份 ODE 建模的检查单，都是从真实事故里提炼的：量纲一致——每个方程左右两边的单位必须对上，导数项多一个时间量纲，漏乘或漏除一个速率常数是最常见笔误；初值完整——几阶方程组就要几个初值，缺一个求解器直接报错；参数取值有出处——速率常数来自文献、实验拟合还是现场估计，要标注，敏感性分析优先扫没把握的那几个；稳态健全性检查——把导数置零解出的平衡态是否符合物理直觉（温度不该出现负绝对温度、数量不该超过容量），一分钟的手算能拦住大量低级错误。这四条检查通过后，再把方程交给求解器，返工率会显著下降。

本节要点回顾

• 建模三步：选状态量、写变化率依赖关系、校核量纲与初值；
• 欧拉一阶、龙格-库塔四阶：误差阶决定"步长减半、误差降多少"；
• 自适应步长是生产默认，等间隔输出靠稠密插值而非拆区间循环；
• 刚性识别：快慢时间尺度共存时显式法步长被压死，换隐式 BDF；
• 解析解能得时先求它：数值结果的校核基准比什么都值钱。