一文搞懂Simulink的定步长求解器

Simulink 提供了定步长（fixed-step）、变步长（variable-step）、离散（discrete）这几类求解器。

默认情况下，Simulink会根据系统的方程类型、状态、连续性等因素自动选择合适的求解器。

对于一个系统而言，我们在选择求解器时一般会考虑以下几点因素：

• 系统动态特性
• 解稳定性
• 计算速度
• 求解器稳健性

在实际应用过程中，由于需求的不同，往往会更加侧重上述几点因素中的某一个，因此自动选择的求解器有时不一定是最合适的。

了解各类求解器的计算原理和优劣，有助于我们根据实际问题灵活选择，达到最好的计算效果。

本文将介绍Simulink中的定步长求解器。

定步长求解器的原理是按固定时间间隔从仿真开始到仿真结束的时间段内解算模型，间隔的大小称为步长。

定步长求解器属于连续型求解器，运用数值积分方法来计算模块定义的连续状态。

定步长求解器类型

定步长求解器的原理是按固定时间间隔从仿真开始到仿真结束的时间段内解算模型，间隔的大小称为步长。

定步长求解器属于连续型求解器，运用数值积分方法来计算模块定义的连续状态。

定步长求解器一共包含8种：

1. ode1（Euler）
2. ode1be (Backward Euler)
3. ode2(Heun)
4. ode3（Bogacki-Shampine）
5. ode4（Runge-Kutta）
6. ode5（Dormand-Prince）
7. ode8（Dormand-Prince）
8. ode14x (extrapolation)

ode1

ode1 求解器使用欧拉积分方法，采用当前状态值和状态导数的显函数来计算下一个时间步的模型状态。

欧拉法是一种一阶显式数值方法，适用于初值问题：

$$\frac{dy}{dt}=f(t,y),y(t_0)=y_0$$

欧拉法的公式为：

$$y_{n+1}=y_n+h\cdot f(t_n,y_n)$$

其中：

• h 是固定步长，
• $t_{n+1}=t_n+h$
• $f(t_n,y_n)$是当前点的导数（斜率）

此求解器比高阶求解器需要更少的计算，但提供的准确性相对较低。

优点：

1. 计算简单：欧拉法是 ODE 数值求解中最基本的算法，易于理解和实现。
2. 计算效率高：每步只需评估一次 ( f(t, y) )，计算成本较低，特别适合大步长场景。
3. 固定步长：作为固定步长求解器，ode1 适合需要恒定时间步长的应用，例如实时仿真或硬件同步。

缺点：

1. 精度较低：欧拉法仅为一阶精度，误差随步长 ( h ) 线性减少（即误差比例于 ( h )，而非 $h^2$或更高）。这意味着为了获得高精度，需要使用非常小的步长，这会增加计算量。
2. 稳定性问题：对刚性（stiff ）问题，欧拉法可能不稳定。
3. 适用范围有限：对于解中有快速变化或振荡的系统，ode1 可能无法准确捕获这些特征，除非使用非常小的步长。

ode1 求解器适合以下问题：

1. 非刚性问题：解变化缓慢，系统动态不剧烈的场景，例如简单的机械系统或线性系统。
2. 简单系统：计算效率比精度更重要的场景，例如初步模型验证或快速原型设计。
3. 实时仿真：需要固定步长以与外部硬件或实时系统同步的应用，例如硬件在环（HIL）测试。

ode1be

ode1be采用后向欧拉（Backward Euler）方法，这是一种隐式数值方法， 其核心思想是通过新时间步的导数来估算解。

其公式为：

$$y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$$

其中：

• ( h ) 是固定的步长，
• $t_{n+1}=t_n+h$
• $f(t_{n+1},y_{n+1})$是新时间步的导数。

由于 $y_{n+1}$出现在右边的函数中，这是一个隐式方程。对于线性问题，可以直接解出 $y_{n+1}$，例如：

$$\begin{array}{c}{y}_{n+1}+hk{y}_{n+1}=y_n\\ y_{n+1}(1+hk)=y_n\\ y_{n+1}=\frac{y_n}{1+hk}\end{array}$$

但对于非线性问题，通常需要使用迭代方法如 Newton 方法求解。

优点:

1. 稳定性：Backward Euler 对线性问题具有无条件稳定性，意味着无论步长 ( h ) 多大，方法都不会发散。
2. 适合刚性问题：由于其稳定性，Backward Euler 允许使用较大的步长来处理快速变化的系统，而不会像显式方法那样需要极小的步长，从而减少计算成本。

缺点:

1. 计算成本：隐式方法需要求解方程，每步可能需要多次迭代，计算成本较高。尤其在非线性问题中，求解过程可能收敛缓慢或需要复杂的数值方法。
2. 实现复杂：相比Euler法的直接计算，Backward Euler 需要额外的求解步骤，可能增加编程和调试的难度。
3. 精度限制：Backward Euler 是一阶精度，误差随步长 ( h ) 线性减少（比例于 ( h )），在非刚性问题中可能不如更高阶方法准确。

与Eular法相比：

例如，考虑 ODE $\frac{dy}{dt}=-y$，初始条件 $y(0)=1$，精确解为 $y(t)=e^{-t}$。使用步长 $h=1$​h = 1：

• Euler：$y_1=1+1\cdot (-1)=0$，与精确值 0.3679 差异较大。
• Backward Euler：$y_1=1+1\cdot (-y_1)$，解得 $y_1=0.5$，结果更接近精确值。

这表明，即使步长较大，Backward Euler的结果更准确。

ode2/3/4

ode2/3/4是固定步长的连续显式求解器，都基于Runge-Kutta 方法，这里放在一起介绍。

ode2采用的Heun 方法（也称为改进的欧拉方法）是一种二阶 Runge-Kutta 方法。

Heun 方法的公式为：

$$y_{n+1}=y_n+\frac{h}{2}(k_1+k_2)$$

其中：

$k_1=f(t_n,y_n$)

$k_2=f(t_n+h,y_n+h{k}_1)$

Heun方法通过取区间开始时的斜率 $k_1$和结束时的斜率 $k_2$ 的平均值，计算下一时间步的值，提供了比一阶欧拉法更高的精度，局部截断误差为 $O(h^3)$，全局误差为 $O(h^2)$。

ode3采用的Bogacki-Shampine 方法是一种三阶 Runge-Kutta 方法，每步约需三次函数求解。其公式为：

$$y_{n+1}=y_n+\frac{2}{9}h{k}_1+\frac{1}{3}h{k}_2+\frac{4}{9}h{k}_3$$

其中：

$k_1=f(t_n,y_n)$

$k_2=f(t_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)$

$k_3=f(t_n+\frac{3}{4}h,y_n+\frac{3}{4}h{k}_2)$

局部截断误差为 $O(h^4)$，全局误差为 $O(h^3)$。

Bogacki-Shampine方法具有 FSAL（First Same As Last）属性，即方法中最后一个阶段的函数值可直接用作下一步的第一个阶段值，从而可以减少计算量，提升效率。

ode4采用经典的四阶 Runge-Kutta 方法（RK4），每步需四次函数求解，其公式为：

$$y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$$

其中：

$k_1=f(t_n,y_n)$

$k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1)$

$k_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2)$

$k_4=f(t_n+h,y_n+h{k}_3)$

局部截断误差为 $O(h^5)$，全局误差为 $O(h^4)$。

ode2/3/4的优点：

1. 更高精度：相比一阶的 ode1，ode2/3/4在相同步长下结果更准确。
2. 总体效率更高：虽然比 ode1 复杂（需要多次函数求解），但其高精度允许使用较大的步长，从而在总体仿真时间上可能更高效。
3. 适合非刚性问题：作为显式方法，ode2 适用于系统动态变化不剧烈的非刚性问题，如简单的机械系统或线性系统。

ode2/3/4的缺点：

1. 稳定性问题：作为显式方法，ode2/3/4对刚性问题可能不稳定。
2. 计算成本较高：相比ode1，ode2/3/4 每步需要多次函数求解（计算 $k_n$），计算成本增加，特别是在高维系统或复杂函数下。
3. 无误差控制：作为固定步长求解器，ode2/3/4 无法根据局部误差动态调整步长，可能导致精度不足（步长过大）或计算浪费（步长过小）。

ode2/3/4 求解器适合以下问题：

1. 非刚性问题：解变化缓慢，系统动态不剧烈的场景，例如简单的机械系统或线性系统。
2. 实时仿真：需要固定步长以与外部硬件或实时系统同步的应用，例如硬件在环（HIL）测试。固定步长确保计算时间可预测，适合实时性要求高的场景。
3. 精度与效率的平衡：当要求比ode1更高的精度时，可以根据计算资源灵活选择ode2/3/4，达到精度与效率的平衡。

ode5/8

ode5/8是定步长连续显式求解器，都基于 Dormand-Prince 方法，这是一种高阶的Runge-Kutta方法。

Dormand-Prince方法通过在每步内计算多个阶段的导数，并以特定权重组合这些值，近似泰勒展开的前几阶项。这种方法在步长区间内多次采样导数，从而显著提高精度。

这里以ode5使用的五阶Dormand-Prince为例，其公式为：

$$y_{k+1}=y_k+\frac{35}{384}{k}_1+\frac{500}{1113}{k}_3+\frac{125}{192}{k}_4-\frac{2187}{6784}{k}_5+\frac{11}{84}{k}_6.$$

其中：

$$\begin{array}{rl}{k}_1& =hf\left(t_k,y_k\right),\\ k_2& =hf\left(t_k+\frac{1}{5}h,y_k+\frac{1}{5}{k}_1\right),\\ k_3& =hf\left(t_k+\frac{3}{10}h,y_k+\frac{3}{40}{k}_1+\frac{9}{40}{k}_2\right),\\ k_4& =hf\left(t_k+\frac{4}{5}h,y_k+\frac{44}{45}{k}_1-\frac{56}{15}{k}_2+\frac{32}{9}{k}_3\right),\\ k_5& =hf\left(t_k+\frac{8}{9}h,y_k+\frac{19372}{6561}{k}_1-\frac{25360}{2187}{k}_2+\frac{64448}{6561}{k}_3-\frac{212}{729}{k}_4\right),\\ k_6& =hf\left(t_k+h,y_k+\frac{9017}{3168}{k}_1-\frac{355}{33}{k}_2-\frac{46732}{5247}{k}_3+\frac{49}{176}{k}_4-\frac{5103}{18656}{k}_5\right),\\ k_7& =hf\left(t_k+h,y_k+\frac{35}{384}{k}_1+\frac{500}{1113}{k}_3+\frac{125}{192}{k}_4-\frac{2187}{6784}{k}_5+\frac{11}{84}{k}_6\right).\end{array}$$

可以看到五阶Dormand-Prince方法有七个阶段，但由于First Same As Last (FSAL) 属性，每步只需六次函数计算，最后阶段的计算（$k_7$）与下一步的 $k_1$相同，可以减少函数评估次数。

以此类推，八阶Dormand-Prince方法需要至少12个阶段，但由于FSAL性质，每步实际只需计算11次函数计算。

ode5/8的优点：

1. 高精度：作为高阶求解器，ode8 的局部截断误差为 $O(h^6)$，全局误差为 $O(h^5)$，而ode8 的局部截断误差为 $O(h^9)$，全局误差为 $O(h^8)$
2. 应用广泛且成熟：Dormand–Prince方法是许多科学计算软件（如 MATLAB、GNU Octave、Python SciPy 等）的默认选择，经过长期应用和大量测试，其数值性能和可靠性已经得到广泛认可。
3. 适合非刚性问题：作为显式方法，Dormand–Prince方法具有较大的稳定区域和良好的收敛性，能够在保证高精度的前提下保持数值稳定。。

ode5/8的缺点

1. 稳定性问题：作为显式方法，ode8 对刚性问题可能不稳定，刚性问题是指解中有快速变化的成分。
2. 无误差控制：作为固定步长求解器，ode5/8 无法根据局部误差动态调整步长，可能导致精度不足（步长过大）或计算浪费（步长过小）。
3. 实现复杂：高阶方法的系数（Butcher 表）较为繁琐，不仅在理论推导上比较复杂，在实际编程实现中也容易出现数值舍入和精度问题，需要仔细处理。
4. 计算成本高：相比低阶求解器，ode5每步需要6次函数计算，ode8 每步需要 13 次函数计算，计算成本较高，特别是在高维系统或复杂函数下。

ode14

ode14x 是一种固定步长隐式求解器，使用Newton方法和外推技术进行积分。其核心思想是将 ODE 视为隐式函数，通过外推提供初始猜测，然后使用 Newton 方法求解隐式方程。

具体来说，ode14x 求解的隐式方程为：

$$X(n+1)=X(n)+h\cdot f(t+h,X(n+1))$$

其中：

• X 是状态向量，
• f 是导数函数（ODE 的右侧函数），
• h 是固定步长，
• t 是当前时间。

由于 $X(n+1)$出现在等式的两侧，这是一个隐式方程，需要迭代方法求解。ode14x 使用 Newton 方法来求解此方程，求解过程分为两个主要步骤：

1. 外推（Extrapolation）:

   使用历史状态和导数信息，计算下一步状态 $X(n+1)$的初始猜测。

   外推阶数由用户设置（1 到 4），高阶外推使用更多历史信息，提供更准确的初始猜测，但计算量更大。

   例如，简单的线性外推可能为 $X(n+1)\approx X(n)+h\cdot f(t,X(n))$，但 ode14x 使用更高级的外推方法，如基于多项式的外推。

2. Newton 方法：

   使用外推得到的初始猜测，迭代求解隐式方程：

   $$X(n+1)-X(n)-h\cdot f(t+h,X(n+1))=0$$

   Newton 方法的迭代公式为：

   $$X_{k+1}=X_k-\frac{F(X_k)}{F^{\prime }(X_k)}$$

   其中：

   $F(X)=X-X(n)-h\cdot f(t+h,X)$

   迭代直到收敛，得到精确的 $X(n+1)$。

这种结合外推和 Newton 方法的策略使得 ode14x 既能提供高精度，又能在固定步长下保持稳定性。

ode14x的优点：

1. 适用于刚性问题：ode14x 作为隐式求解器，适合处理刚性问题（解中有快速变化的成分），显式求解器在这些问题上可能不稳定。
2. 固定步长特性：适合需要固定步长的实时仿真和硬件同步场景，例如硬件在环（HIL）测试。
3. 高精度：通过 Newton 方法和高阶外推，ode14x 提供高精度的解，特别是在固定步长下。

ode14x的缺点：

1. 计算成本高：每步需要多次 Newton 迭代和外推计算，计算量较大，特别是在高维系统或复杂函数下。
2. 无误差控制：作为固定步长求解器，ode14x 无法根据局部误差动态调整步长，可能导致精度不足（步长过大）或计算浪费（步长过小）。
3. 调试复杂：由于隐式方法和迭代过程，调试可能比显式求解器更复杂。

总结

本文介绍了Simulink中的8种定步长求解器，包括6种显式求解器，2种隐式求解器，其差异主要在于精度和计算成本上，如下表所示：

求解器	方法	显式/隐式	精度（阶数）	计算成本
ode1	Euler’s method	显式	1st	低
ode1be	Backward Euler	隐式	1st	高（固定）
ode2	Heun’s method	显式	2nd	中等
ode3	Bogacki-Shampine	显式	3rd	中等至高
ode4	Runge-Kutta 4th order (RK4)	显式	4th	高
ode5	Dormand-Prince 5th order	显式	5th	高
ode8	Dormand-Prince 8th order	显式	8th	极高
ode14x	Newton + extrapolation	隐式	最高达4th	高（迭代）

关于文中涉及的“刚性” 这个概念，可参考数值计算与仿真中的 “刚性” 是什么？一文

---