📚 课程目录第4章数值仿真基础与工具选择

4. 数值仿真基础与工具选择

Lorenz系统是一组三维非线性常微分方程，除少数参数组合外没有解析解，因此数值积分是研究其混沌行为的唯一可行途径。本节围绕三个问题展开：用什么方法积分、如何判断方法是否可靠、以及如何用Python或MATLAB搭建一个可复用的仿真框架。

一、从连续方程到离散迭代

Lorenz系统的标准形式为：

dx/dt = σ(y − x)
dy/dt = x(ρ − z) − y
dz/dt = xy − βz

经典参数取 σ = 10、ρ = 28、β = 8/3，此时系统呈混沌状态。数值积分的基本思路是：已知 t_n 时刻的状态 u_n = (x_n, y_n, z_n)，用某种规则估计区间 [t_n, t_n+h] 内的平均斜率，从而推进到 u_n+1。步长 h 的选择直接决定精度与稳定性。

二、欧拉法与龙格-库塔法

显式欧拉法只取区间起点的斜率：u_n+1 = u_n + h·f(u_n)。它实现简单，但只有一阶精度，局部截断误差为 O(h²)。对Lorenz系统而言，欧拉法的实际问题是：当 h 稍大（例如 h = 0.01 以上）时，数值轨迹会迅速偏离真实轨道，甚至出现能量虚假增长导致发散。因此欧拉法只适合教学演示，不适合定量研究。

四阶龙格-库塔法（RK4）在每个步长内计算四次斜率并加权平均：

k1 = f(u_n)
k2 = f(u_n + h/2 · k1)
k3 = f(u_n + h/2 · k2)
k4 = f(u_n + h · k3)
u_{n+1} = u_n + h/6 · (k1 + 2k2 + 2k3 + k4)

RK4 具有四阶精度，局部误差 O(h⁵)，在相同步长下精度远高于欧拉法，是Lorenz系统仿真的常用默认选择。代价是每步四次函数求值，计算量约为欧拉法的四倍，但对三维系统而言这一开销可以忽略。

三、精度与稳定性的比较

需要区分两个概念：精度指单步或全局误差随步长缩小的阶数；稳定性指误差是否会被放大。Lorenz系统本身对初值极端敏感（正Lyapunov指数约0.9），任何数值方法都无法长期保持轨迹逐点一致，因此实践中关注的是：短时间内的轨迹精度，以及长时间统计量（如吸引子形状、均值、功率谱）是否收敛。

一个实用的验证方法是步长收敛测试：用 h、h/2、h/4 分别积分同一初值到固定时刻 T，比较终点状态差。若差值随步长减半而按 2⁴ 比例缩小，说明RK4进入四阶收敛区；若差值不缩小或反而增大，说明步长仍偏大或已进入混沌放大区。对Lorenz系统，RK4 常用 h = 0.001～0.01，配合 T = 20～50 的积分区间即可得到形状正确的吸引子。

判断步长是否合适的经验做法：先用 h = 0.01 与 h = 0.005 各积分到 T = 20，若两条轨迹在 t < 10 内肉眼重合、之后才因混沌分离，说明步长已足够解析局部动力学；若在 t < 5 就明显分叉，应减小步长。

四、工具选择：Python 与 MATLAB

两者都能胜任Lorenz仿真，差异主要在生态与使用习惯。Python 的优势是免费、库丰富：scipy.integrate.solve_ivp 内置 RK45（自适应步长的Runge-Kutta-Fehlberg）、RK23、DOP853 等方法，配合 NumPy 做向量化、Matplotlib 做三维绘图，适合批量参数扫描与后续机器学习衔接。MATLAB 的优势是 ode45 开箱即用、矩阵语法直观、Simulink 便于搭建框图式模型，在控制系统教学与工程验证中更常见。

选择建议：若研究重点是参数扫描、分岔图、Lyapunov指数计算等需要大量重复积分的任务，Python 的脚本化与并行能力更顺手；若侧重课堂演示、与控制系统工具箱联动，MATLAB 更直接。两者算法本质相同，结论可互相验证。

五、搭建基础仿真框架

一个可复用的框架应包含四部分：方程定义、积分器调用、结果存储、可视化与验证。以下为 Python 示例（参数为经典值，结果为示意性输出，需实际运行确认）：

import numpy as np
from scipy.integrate import solve_ivp

def lorenz(t, u, sigma=10.0, rho=28.0, beta=8/3):
    x, y, z = u
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

u0 = [1.0, 1.0, 1.0]
t_span = (0, 40)
t_eval = np.linspace(0, 40, 4001)

sol = solve_ivp(lorenz, t_span, u0, method='RK45',
                t_eval=t_eval, rtol=1e-8, atol=1e-10)
print(sol.y.shape)   # (3, 4001)
print(sol.y[:, -1])  # 终点状态

关键设置说明：rtol 与 atol 控制自适应步长的误差容限，混沌系统建议收紧到 1e-8 量级；t_eval 指定输出采样点，避免自适应步长导致输出不均匀。若改用固定步长RK4，可自行实现上面的四级公式，便于与自适应方法对照。

MATLAB 对应写法为：

f = @(t,u) [10*(u(2)-u(1));
            u(1)*(28-u(3))-u(2);
            u(1)*u(2)-8/3*u(3)];
[t, U] = ode45(f, [0 40], [1;1;1]);
plot3(U(:,1), U(:,2), U(:,3));

框架搭建完成后，应做两项基本检查：其一，用两种方法（如 RK4 与 ode45）积分同一初值，比较吸引子形态是否一致；其二，改变初值 1e-6 量级，观察轨迹分离时间是否符合混沌预期。这两步能有效排除步长过大、参数写错等常见错误。

最后需要强调：混沌系统的数值仿真不存在“绝对正确”的长期轨迹，任何方法都会在有限时间后与参考解分离。因此评价仿真质量的标准不是逐点吻合，而是吸引子几何形态、统计分布与参数响应趋势是否稳定。把步长收敛测试和双方法交叉验证纳入固定流程，比单纯追求更小步长更有意义。