第四章：连续系统仿真——微分方程建模与数值积分方法

各位工程师朋友，欢迎来到连续系统仿真这一章。说实话，连续系统仿真在工程领域太常见了——从电路中的RC充放电，到机械系统的弹簧阻尼，再到飞行器的姿态控制，本质上都是在解微分方程。我个人觉得，搞懂这一章，你就能看懂大半的物理仿真问题。

4.1 微分方程建模——把物理世界翻译成数学语言

先聊聊建模。微分方程建模，说白了就是把物理规律写成数学等式。比如牛顿第二定律 F=ma，写成微分形式就是：

m * d²x/dt² = F(t) - b * dx/dt - k * x

这里 x 是位移，dx/dt 是速度，d²x/dt² 是加速度。b 是阻尼系数，k 是弹簧刚度。你看，一个机械振动系统就这么出来了。

建模的核心步骤，我总结为三步：

1. 确定状态变量——系统里哪些量会随时间变化？比如位置、速度、温度、电压。
2. 写出守恒定律——能量守恒、动量守恒、质量守恒，这是物理基础。
3. 整理成标准形式——把高阶微分方程降阶为一阶方程组，方便计算机处理。

4.2 数值积分方法——计算机怎么解微分方程？

微分方程写好了，但解析解往往求不出来。怎么办？用数值方法。计算机不会解方程，但它会做加法。数值积分就是把这个思想用到微分方程求解上。

考虑一阶微分方程：

dy/dt = f(t, y),  y(t₀) = y₀

我们要从 t₀ 开始，一步步算出 y(t₁), y(t₂), ... 直到终点。步长 h = t_{n+1} - t_n 决定了精度和速度。

4.2.1 欧拉法——最简单，但别小看它

欧拉法的思路很直白：用当前点的斜率去预测下一个点的值。

y_{n+1} = y_n + h * f(t_n, y_n)

代码实现也就几行：

def euler(f, y0, t0, tf, h):
    t = t0
    y = y0
    while t < tf:
        y = y + h * f(t, y)
        t = t + h
    return y

嗯，这里要注意——欧拉法虽然简单，但精度只有 O(h)。什么意思？步长减半，误差只减半。对于刚性系统或者长时间仿真，误差会累积到不可接受。

4.2.2 龙格-库塔法——工程界的标配

龙格-库塔法（RK）是欧拉法的升级版。它不是在当前点取一个斜率，而是在区间内取多个点的斜率，然后加权平均。最常用的是四阶龙格-库塔法（RK4）。

RK4 的计算公式：

k₁ = f(t_n, y_n)
k₂ = f(t_n + h/2, y_n + h*k₁/2)
k₃ = f(t_n + h/2, y_n + h*k₂/2)
k₄ = f(t_n + h, y_n + h*k₃)
y_{n+1} = y_n + (h/6) * (k₁ + 2k₂ + 2k₃ + k₄)

代码实现：

def rk4(f, y0, t0, tf, h):
    t = t0
    y = y0
    while t < tf:
        k1 = f(t, y)
        k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2)
        k4 = f(t + h, y + h*k3)
        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
    return y

RK4 的精度是 O(h⁴)。步长减半，误差减到原来的 1/16。这就是为什么它成为工程仿真的首选。

4.3 系统动力学仿真——把方法用到真实系统

理论讲完了，咱们看个实际例子。考虑一个简单的 RLC 电路，或者一个质量-弹簧-阻尼系统。它们的微分方程形式是一样的：

d²x/dt² + 2ζω₀ dx/dt + ω₀² x = F(t)/m

其中 ω₀ 是自然频率，ζ 是阻尼比。降阶为一阶方程组：

dx₁/dt = x₂
dx₂/dt = -2ζω₀ x₂ - ω₀² x₁ + F(t)/m

用 RK4 仿真这个系统：

import numpy as np

def system(t, y, zeta, omega0, m):
    x1, x2 = y
    F = 1.0  # 阶跃输入
    dx1 = x2
    dx2 = -2*zeta*omega0*x2 - omega0**2*x1 + F/m
    return np.array([dx1, dx2])

# 参数设置
zeta = 0.2
omega0 = 2.0 * np.pi
m = 1.0
y0 = np.array([0.0, 0.0])
t0, tf, h = 0.0, 5.0, 0.01

# 仿真循环
t = t0
y = y0
while t < tf:
    y = rk4(lambda t, y: system(t, y, zeta, omega0, m), y, t, t+h, h)
    t = t + h
    print(f"t={t:.3f}, x={y[0]:.5f}, v={y[1]:.5f}")

你想想看，这个简单的代码背后，其实涵盖了连续系统仿真的全部核心：建模、降阶、数值积分、步长控制。我在做汽车悬架仿真时，用的就是这套框架，只是把参数换成了更复杂的非线性阻尼模型。

4.4 知识体系总览

为了让你更直观地理解本章的知识结构，我画了一张图：

4.5 方法对比与选型建议

不同方法各有优劣，我整理了一个对比表：

方法	精度阶数	每步计算量	稳定性	推荐场景
欧拉法	O(h)	1次函数求值	较差	快速原型、教学演示
改进欧拉法	O(h²)	2次函数求值	一般	精度要求不高的仿真
RK4	O(h⁴)	4次函数求值	良好	工程仿真首选
自适应RK45	O(h⁴)~O(h⁵)	6次函数求值	优秀	变刚度系统、长时间仿真

好了，连续系统仿真的核心内容就这些。微分方程建模是基础，数值积分方法是工具，系统动力学仿真是应用。三者缺一不可。你在实际项目中遇到什么问题，欢迎随时交流。

---