Part I · 第 2 章 微积分 · 2.5

常微分方程

微分方程不告诉你"状态是多少"，而是告诉你在每个状态下它会怎样变化。电机转速、关节角度、电池电量、机器人位置，凡是随时间演化的东西都由它描述。解微分方程就像顺着风向图走：从起点出发，每一刻都朝箭头指的方向走一小步。

约 30 分钟 前置：2.1 导数 · 2.3 积分 · 1.11 特征值（推荐） 5 个交互

    1为什么机器人需要它

    机器人是一个动力系统：给它一个输入（电压、力矩），它的状态按物理规律随时间演化。写成公式就是 \(\dot{\mathbf x}=\mathbf f(\mathbf x,\mathbf u)\)：状态的变化率由当前状态和输入决定。

    • 电机的转速、关节的角度、足式机器人的质心位置，都服从微分方程。
    • 控制器的设计，就是挑选输入 \(\mathbf u\)，让这个方程的解沿着我们想要的方向走，而且稳定。
    • 仿真器（MuJoCo、Isaac、PyBullet）每一步都在数值求解一个微分方程；强化学习里用到的"环境转移"，在连续的世界里就是它的离散化。

    这一节的目标：读懂 \(\dot x=f(x)\) 的几何意义，会判断一个线性系统稳不稳，理解为什么步长和方法的选择会让仿真失效。

    2直觉：顺着箭头走

    \(\dot x=f(t,x)\) 的意思是：在平面上的每一点 \((t,x)\)，都规定了一个斜率。把这些斜率画成一片小线段，就是斜率场（方向场）。解 \(x(t)\) 就是一条处处与小线段相切的曲线：它像一片叶子顺着水流漂。换一个起点，就是另一条解，而同一点只会有一条解（在合理的条件下），所以解的曲线不会相交。

    在下面的图里，拖动黄点（初始条件），点击空白处再放几条解。换几个方程看看。

    FIG 1斜率场与解
    小线段是各点的斜率，琥珀色曲线是可拖动初始条件的主解（向前向后都算），青色是点击放下的其他解。注意：曲线一旦"吸附"到某个水平的位置，说明那里是平衡点（斜率为 0）；稳定的平衡点把附近的解都吸引过来。
    • 选"逻辑斯蒂"，在 \(x>0\) 处放下至少 3 条不同起点的解（3 的上方和下方都试试）：它们都靠向同一条水平线（稳定平衡点 \(x=3\)）。
    • 选"\(\dot x=0.25x^2\)"，把起点放到 \(x_0>1\)：解在有限时间内爆炸，曲线冲出画面。
    核心洞见

    微分方程把"随时间怎样变"写成"当前状态的函数"。解不是公式，而是一条被箭头牵着走的轨迹。有了这个图景，"稳定"就是"所有附近的轨迹都靠向这里"，"数值积分"就是"沿着箭头一小步一小步走"。

    3动手玩：一阶、二阶、数值方法

    一阶系统：电机转速

    一个带摩擦的电机，转速 \(\omega\) 满足 \(\tau\,\dot\omega=-\omega+K u\)：\(\tau\) 是时间常数，\(K\) 是稳态增益。输入阶跃 \(u=1\)，解是 \(\omega(t)=K\,(1-e^{-t/\tau})\)：每过一个 \(\tau\)，离终值的距离缩小到原来的 \(1/e\approx37\%\)；过 \(\tau\) 时到达终值的 63%，过 \(4\tau\) 时到 98%。输入是正弦时，系统像一个低通滤波器：高频被削弱，还会滞后。

    FIG 2一阶响应
    琥珀色是输入 \(u(t)\)，青色是转速 \(\omega(t)\)（用 RK4 数值求解）。竖虚线标出 \(t=\tau\)，横虚线标出终值的 63%。
    • 阶跃输入下，调时间常数 \(\tau=0.5\) s：转速在 \(t=0.5\) s 时到达终值的 63%。
    • 切到"正弦"，调频率使 \(\omega_{in}\tau\approx1\)：输出幅值降到约 71%，滞后 45°。

    二阶系统：质量–弹簧–阻尼器，与极点

    关节、悬架、PD 控制的机械臂，常见形式是 \(\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=\omega_n^2 r\)。\(\omega_n\) 是固有频率，\(\zeta\) 是阻尼比。把它写成 \(\dot{\mathbf x}=A\mathbf x\)（状态 \(\mathbf x=(\theta,\dot\theta)\)），\(A\) 的特征值（又叫极点）决定一切：解里含 \(e^{\lambda t}\)，实部 \(\mathrm{Re}\,\lambda<0\) 才会衰减，虚部决定振荡。右图是复平面，拖动紫色的极点看看响应怎么变；极点越靠左越快，越靠近虚轴越慢，跨过虚轴就发散。

    FIG 3阻尼比与极点
    阶跃响应 θ(t)
    极点（复平面）
    左：对单位阶跃输入的响应。右：极点位置，紫色点可拖动（\(\zeta<1\) 时是一对共轭极点；\(\zeta\ge1\) 时两个极点落在实轴上）。灰色半圆是 \(|\lambda|=\omega_n\)。
    • 找到临界阻尼 \(\zeta=1\)：最快且不超调的响应。
    • 让超调量约为 16%（\(\zeta\approx0.5\)，工程里常用的折中）。
    • 把阻尼调成负数（拖动极点到虚轴右侧）：响应越摆越大，系统不稳定。

    数值方法：步长、方法与"仿真爆炸"

    一个理想的无阻尼振子 \(\ddot x=-x\)，其能量 \(E=\tfrac12(x^2+v^2)\) 应该永远不变，相图里的轨迹是一个圆。用不同的数值方法、不同的步长 \(h\) 积分，看谁能保住这个圆。

    FIG 4谁能保住能量？
    相图 (x, v)
    log₁₀ (E / E₀) 随时间
    灰色圆是精确解。显式欧拉每步把能量乘上 \((1+h^2)\)，轨迹向外螺旋；半隐式（辛）欧拉能量在一个小范围内振荡，永不漂走；RK4 精度高，但在大步长下能量也会缓慢衰减。
    • 选显式欧拉，步长 \(h\ge0.2\)：20 秒内能量涨到初始的 2 倍以上。
    • 选半隐式欧拉，同样的步长：能量保持在 ±25% 以内，没有漂走。

    4再看数学

    提示：悬停在公式带底色的符号上，FIG 3 中对应的部分会高亮。

    标准形式与状态空间

    任何高阶方程都可以写成一阶方程组。例如 \(\ddot\theta=g(\theta,\dot\theta,u)\)：令 \(\mathbf x=(\theta,\dot\theta)\)，则 \(\dot{\mathbf x}=(\dot\theta,\ g(\theta,\dot\theta,u))\)。状态 \(\mathbf x\) 是"知道它就能预测未来"的最少信息。

    $$ \dot{\mathbf x}=\mathbf f(\mathbf x,\mathbf u,t),\qquad \mathbf x(t_0)=\mathbf x_0. $$

    只要 \(\mathbf f\) 光滑（局部 Lipschitz），给定初值，解存在并且唯一：所以解曲线不相交。但解可能在有限时间内爆炸（FIG 1 的 \(\dot x=0.25x^2\)）。

    线性系统与极点

    $$ \dot{\mathbf x}=A\mathbf x\ \Longrightarrow\ \mathbf x(t)=e^{At}\mathbf x_0=\sum_i c_i\,e^{\hl{lam}{\lambda_i}t}\,\mathbf v_i,\qquad A\mathbf v_i=\lambda_i\mathbf v_i. $$

    沿每个特征向量 \(\mathbf v_i\)，系统退化成一个一阶方程 \(\dot c=\lambda c\)，解是 \(e^{\lambda t}\)（特征值的概念见 1.11 节）。所以：\(\mathrm{Re}\,\lambda_i<0\) 对所有 \(i\) 成立，则系统渐近稳定；只要有一个 \(\mathrm{Re}\,\lambda_i>0\)，就不稳定；虚部给出振荡频率。二阶系统 \(\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=0\) 的极点是

    $$ \lambda_{1,2}=\hl{sigma}{-\zeta\omega_n}\;\pm\;\hl{wd}{\,j\,\omega_n\sqrt{1-\zeta^2}}\qquad(\zeta<1),\qquad M_p=e^{-\pi\zeta/\sqrt{1-\zeta^2}}. $$

    其中 \(M_p\) 是阶跃响应的超调量（峰值超出终值的比例）：\(\zeta=0.5\) 时约 16%，\(\zeta=0.7\) 时约 5%，\(\zeta\ge1\) 时为 0。

    一阶系统 \(\tau\dot\omega=-\omega+Ku\) 的极点是 \(-1/\tau\)：时间常数越小，极点越靠左，响应越快。对非线性系统，在平衡点处用泰勒展开（2.4 节）得到线性化 \(\dot{\delta\mathbf x}=J\,\delta\mathbf x\)，\(J\) 是雅可比，它的特征值给出这个平衡点的局部稳定性。

    数值求解

    方法一步阶
    显式欧拉\(\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f(\mathbf x_k)\)1
    中点法（RK2）\(\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f\big(\mathbf x_k+\tfrac h2\mathbf f(\mathbf x_k)\big)\)2
    RK4\(\mathbf x_{k+1}=\mathbf x_k+\tfrac h6(\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4)\)4
    半隐式欧拉（辛）\(v_{k+1}=v_k+h\,a(x_k),\ x_{k+1}=x_k+h\,v_{k+1}\)1

    它们都是 2.3 节积分规则的推广（欧拉 = 左矩形，中点法 = 中点规则，RK4 ≈ Simpson）。还有一个 2.3 节没有的问题：稳定性。对 \(\dot x=\lambda x\)（\(\lambda<0\)），显式欧拉要求 \(h<2/|\lambda|\)，否则数值解自己会发散，哪怕真解在衰减。系统里同时存在很快和很慢的模态（"刚性"）时，只能用极小的步长，或改用隐式方法。

    5机器人中的它：PD 控制的关节

    一根受重力作用的关节：\(J\ddot\theta+b\dot\theta+mg\ell\sin\theta=\tau\)。控制器输出力矩 \(\tau=k_p(\theta^\star-\theta)-k_d\dot\theta\)（PD 控制）。线性化后（并补偿重力时）就是 FIG 3 的二阶系统：\(\omega_n=\sqrt{k_p/J}\)，\(\zeta=\dfrac{k_d+b}{2\sqrt{k_pJ}}\)。\(k_p\) 像弹簧，决定"多快、多硬"；\(k_d\) 像阻尼，抑制振荡。下面让它去追一个 60° 的目标，看重力怎样让它停在错误的位置。

    FIG 5PD 控制的关节
    参数：\(J=0.5\) kg·m²，\(b=0.1\) N·m·s/rad，\(mg\ell=4\) N·m。曲线是完整非线性方程的数值解（RK4），灰色带是 ±2% 的目标范围。
    • 不补偿重力，\(k_p\le3\)：关节停在离目标超过 10° 的地方（稳态误差）。再打开重力补偿，误差消失。
    • 打开重力补偿，调 \(k_p\)、\(k_d\)，使超调 小于 5%，且 1.5 秒内进入 ±2% 范围。
    真实系统里

    真实控制器是在离散时间里运行的：每个周期（如 1 kHz）读传感器、算力矩、写给电机。这相当于对上面的连续方程做了一次"零阶保持"的离散化；\(k_p\) 太大、周期太长，闭环极点会越过离散系统的稳定边界，机械臂会开始高频抖动。所以控制器的采样周期要比时间常数小得多，这和上面欧拉法的步长限制是同一件事。后续 Part IV 会系统地讲 PID、LQR 和 MPC。

    6常见误区

    误区 1：步长越大越快，所以用大步长

    步长必须小于系统最快模态的时间尺度。显式欧拉在 \(h>2/|\lambda|\) 时发散，而且常常不是"算不准"，是"炸掉"：数值变成 NaN。仿真里突然出现的"机器人飞走了"，很多就是这个。

    误区 2：欧拉法对振荡系统也够用

    显式欧拉每一步都往能量里"注入"一点，纯振荡的系统会越摆越大。刚体动力学里常用半隐式欧拉或更高阶方法。

    误区 3：只看稳不稳，不看快不快、抖不抖

    所有极点都在左半平面只保证最终会收敛。靠近虚轴的极点收敛得极慢；虚部很大的极点收敛时伴随强烈的振荡。设计时要同时看"实部（快慢）"和"阻尼比（抖不抖）"。

    误区 4：线性化的稳定性适用于一切

    线性化只在平衡点附近有效（2.4 节）。FIG 5 里不补偿重力时，系统停在的并不是目标，而是另一个平衡点。非线性系统可以有多个平衡点，每个都有自己的局部稳定性。

    7小练习

    1. \(\dot x=-2x\)，\(x(0)=3\)。\(x(1)\) 约是多少？
    解是 \(x(t)=3e^{-2t}\)，\(x(1)=3e^{-2}\approx0.406\)。
    2. 一阶系统时间常数 \(\tau=0.2\) s，对单位阶跃输入，\(t=0.2\) s 时输出大约到终值的多少？
    \(1-e^{-1}\approx0.632\)。到 \(2\tau\) 时是 86%，到 \(4\tau\) 时是 98%。
    3. 线性系统的极点是 \(-1\pm2j\)。它的行为是？
    实部 −1 < 0：稳定，衰减的包络是 \(e^{-t}\)；虚部 ±2：振荡角频率是 2 rad/s。
    4. 二阶系统阻尼比 \(\zeta=1\)（临界阻尼），阶跃响应是？
    两个极点重合在 \(-\omega_n\)。\(\zeta<1\) 会超调，\(\zeta>1\) 的慢极点更靠近原点，响应更慢。
    5. 用显式欧拉模拟一个无阻尼振子，时间一长会看到什么？
    显式欧拉每步把能量乘 \(1+h^2\)，指数增长。半隐式欧拉则没有这个漂移。

    8延伸阅读

    • 3Blue1Brown, Differential Equations — 第 1 章、第 2 章：相空间与方向场
    • S. Strogatz, Nonlinear Dynamics and Chaos — 一维、二维流与稳定性
    • K. Åström & R. Murray, Feedback Systems — 第 3、5 章：动态行为与线性系统
    • E. Hairer, S. Nørsett & G. Wanner, Solving Ordinary Differential Equations I — 数值方法的经典教材
    • K. Lynch & F. Park, Modern Robotics, Ch. 11 — 机器人控制：二阶系统与 PD 控制