《飞控数学与控制 · 新手入门连载》第 5 期｜数值离散化：欧拉法与 RK4——让微分方程跑进单片机

坐看风云变幻

316人浏览 · 2026-10-14 11:18:54

坐看风云变幻 · 2026-10-14 11:18:54 发布

《飞控数学与控制 · 新手入门连载》第 5 期｜数值离散化：欧拉法与 RK4——让微分方程跑进单片机

导读：前四期我们把无人机写成连续微分方程（第 3 期），又知道测量带噪声（第 4 期）。但飞控跑在单片机里，状态只能按固定采样周期一步步更新，不是连续流动。本期讲"连续 → 离散"的核心方法：欧拉法怎么从导数定义直接"抄"出迭代式，RK4 为什么更准，采样率该选多高。这是把第 3 期模型变成可运行代码的关键一步，也是姿态解算、位置积分、滤波实现的通用工具。全程零基础可读，公式逐步推导不跳步。

一、为什么飞控代码里没有"微分方程"？

第 3 期的动力学方程 $x˙=f(x,t)\dot{x} = f(x, t)$ 是连续的：时间 $t$ 连续变化。但飞控代码是这样跑的：

while (true) {
    读取传感器；          // 每 1ms 或 4ms
    更新状态估计；        // 姿态解算、积分
    计算控制量；          // PID
    输出 PWM；            // 驱动电机
    delay(采样周期);
}

状态不是连续演化的，而是每隔 $Δt\Delta t$ 跳一步。把连续方程改写为"每步更新"的迭代式，就是离散化（数值积分）。

直观理解：连续模型是"水流"，代码是"按帧拍照"。每帧间隔就是采样周期 $Δt\Delta t$ 。图 5-1 展示连续曲线与离散采样的关系。

图5-1 连续解与离散采样：真实状态是光滑曲线，飞控只在采样时刻观测/更新，采样点越密（越小）越接近真实曲线

二、欧拉法：从导数定义直接"抄"

2.1 导数定义的启发

第 3 期我们写过导数的定义：

$\dot{x}(t) = \lim_{\Delta t \to 0}\frac{x(t+\Delta t) - x(t)}{\Delta t}$

技巧：极限近似——飞控的 $Δt\Delta t$ 很小但不为零，把极限符号拿掉：

$\dot{x}(t) \approx \frac{x(t+\Delta t) - x(t)}{\Delta t}$

2.2 前向欧拉迭代式

由 $x˙=f(x,t)\dot{x} = f(x, t)$ ，代入上式（技巧：解出 $x(t+Δt)x(t+\Delta t)$ ）：

$\frac{x(t+\Delta t) - x(t)}{\Delta t} = f(x(t), t)$

$x(t+\Delta t) = x(t) + \Delta t \cdot f(x(t), t)$

写成离散记号 $xk=x(tk)t_k = k\Delta t,\ x_k = x(t_k)$ ：

$\boxed{\,x_{k+1} = x_k + \Delta t \cdot f(x_k, t_k)\,}$

这就是前向欧拉法——飞控里最常用、最简单的数值积分。

飞控实例：姿态角积分。设姿态角 $θ\theta$ 与角速度 $ω\omega$ 的关系 $θ˙=ω\dot{\theta} = \omega$ ，欧拉法直接给出：

$\theta_{k+1} = \theta_k + \Delta t \cdot \omega_k$

这正是姿态解算里"角速度积分得角度"的代码原型。

三、欧拉法的误差：为什么不够用

3.1 截断误差的直观来源

欧拉法用区间起点的斜率代表整个区间，如图 5-2 所示。当曲线弯曲时，起点斜率偏离平均斜率，产生误差。

图5-2 欧拉法 vs 真实解：欧拉法用区间起点的斜率推进（红色折线），真实解是弯曲曲线（蓝色）。每步的竖直差即局部截断误差，步长越大误差越大

3.2 误差阶次

把真实解做泰勒展开（技巧：泰勒展开 + 局部截断）：

$x(t+\Delta t) = x(t) + \Delta t\, f + \frac{\Delta t^2}{2}\ddot{x} + O(\Delta t^3)$

欧拉法只保留到一阶项（技巧：丢掉高阶项）：

$x_{k+1} = x_k + \Delta t\, f(x_k, t_k)$

局部截断误差（每步误差）：

$e_{\text{local}} = \frac{\Delta t^2}{2}\ddot{x} = O(\Delta t^2)$

全局误差（技巧：误差累积——N 步累加， $\propto 1/\Delta t$ ）：

$e_{\text{global}} = N \cdot O(\Delta t^2) = O(\Delta t)$

结论：欧拉法是一阶精度——步长减半，误差大致减半。飞控里 4ms 控制周期积分姿态，欧拉法误差通常在可接受范围；但需要高精度时（如高速旋转、长时间积分）就要用高阶方法。

四、RK4：四阶龙格-库塔法

RK4 的思想（技巧：斜率加权平均）：在区间内取 4 个斜率点，加权平均作为"代表斜率"，精度大幅提升。

4.1 四个斜率

设当前状态 $x_k$ 、时间 $t_k$ ，步长 $\Delta t$ ：

$k_1 = f(x_k, t_k)$

$k_2 = f\left(x_k + \frac{h}{2} k_1,\ t_k + \frac{h}{2}\right)$

$k_3 = f\left(x_k + \frac{h}{2} k_2,\ t_k + \frac{h}{2}\right)$

$k_4 = f(x_k + h k_3,\ t_k + h)$

技巧：中点试探—— $k_2$ 、 $k_3$ 用中点斜率试探改进； $k_4$ 用终点试探。

4.2 加权组合

$x_{k+1} = x_k + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$

技巧：Simpson 权重 1:2:2:1——这是四阶精度的关键组合，局部截断误差 $O(h^5)$ 、全局误差 $O(h^4)$ 。

4.3 为什么权重是 1:2:2:1？

把 RK4 与真实解的四阶泰勒展开对比（技巧：泰勒展开系数匹配）：

$\frac{h^2}{2}f' + \frac{h^3}{6}f'' + \frac{h^4}{24}f''' + \cdots$

把 $k1∼k4k_1 \sim k_4$ 各自展开代入组合式，逐项对比系数，恰好匹配到 $h^4$ 项， $h^5$ 项成为首个不匹配项——所以是四阶精度。权重 1:2:2:1 正是这一匹配的结果，不是经验凑数。

飞控应用：高精度位置积分、轨迹预测、SITL 仿真中常用 RK4；控制律计算（如姿态环）为省算力仍用欧拉法。图 5-3 对比欧拉法与 RK4 的精度。

图5-3 欧拉法 vs RK4 数值积分对比：相同步长下，RK4 几乎贴合真实解，欧拉法明显偏移。步长减小时两者误差都缩小，但欧拉法收敛慢（一阶）、RK4 收敛快（四阶）

五、采样率选择：多快才算够？

5.1 Nyquist 准则

对频率为 $fmaxf_{\text{max}}$ 的信号，采样率必须满足（技巧：香农采样定理）：

$f_s \ge 2 f_{\text{max}}$

飞控实践：姿态环带宽典型 10~30 Hz，控制采样率通常取 250~500 Hz（ $Δt=2\Delta t = 2$ ~4 ms），远超 Nyquist 下限，主要受执行器响应（电机）和算力约束。

5.2 采样率与数值精度

采样率提高（ $Δt\Delta t$ 减小）带来三重收益（技巧：误差-步长关系）：

欧拉法全局误差 $O(Δt)O(\Delta t)$ ：采样率翻倍 → 误差减半；
更早发现高频扰动（如桨叶不平衡振动）；
离散化延迟更小，控制环相位裕度更好。

但也付出代价：CPU 占用、传感器噪声放大。工程平衡点通常在"满足带宽需求的 10~20 倍"。

5.3 实践对照表

环节	典型带宽	常用采样率	积分方法
姿态内环（角速度→力矩）	30~100 Hz	250~1000 Hz	欧拉法
姿态外环（角度→角速度）	10~30 Hz	100~250 Hz	欧拉法
位置环	1~5 Hz	50~100 Hz	RK4（可选）
导航解算（SITL/仿真）	—	任意	RK4

六、小结与下期预告

小结：

飞控代码按采样周期逐步更新，连续方程必须离散化；
欧拉法直接来自导数定义的极限近似，简单但只有一阶精度；
RK4 用 1:2:2:1 加权 4 个斜率，达到四阶精度，适合高精度积分；
采样率选择遵循 Nyquist 准则，兼顾精度、算力与噪声。

思考题：

用欧拉法离散 $v˙=a\dot{v} = a$ （ $a$ 常数），写出迭代式。（答案： $vk+1=vk+aΔtv_{k+1} = v_k + a\Delta t$ ）
RK4 权重为什么是 1:2:2:1？与泰勒展开的系数匹配有什么关系？
姿态环带宽 50 Hz，采样率 250 Hz 满足 Nyquist 吗？裕量是多少倍？（答案：满足，5 倍）