数值计算与仿真中的 “刚性” 是什么？

最新推荐文章于 2024-12-16 14:49:48 发布

充电君

最新推荐文章于 2024-12-16 14:49:48 发布

阅读量1.4k

点赞数 25

分类专栏： AI 文章标签：算法基算仿真

本文链接：https://blog.csdn.net/CDB3399/article/details/143133425

版权

AI 专栏收录该内容

12 篇文章

订阅专栏

拉到文末有惊喜

在数值计算分析领域，刚性（Stiffness） 描述的是在求解微分方程时，由于系统中存在显著不同的时间尺度或空间尺度，导致数值解的稳定性和计算效率出现问题的现象。

刚性一般用于描述包含多个时间尺度的常微分方程（ODE），由于存在多个不同的时间尺度，有的部分变化很快，有的部分变化很慢。这会导致传统的数值方法（如显式欧拉法）需要非常小的步长才能保持数值稳定性。

刚性本身并不是一个数学上严格定义的术语，但它在数值计算及模拟中具有重要的应用意义。其可类比于物理学中的刚体概念。就像刚体对外力会产生瞬时反应一样，刚性方程或系统中的某些部分对变化的响应也极其迅速。

在数值计算与模拟中，刚性的概念一般用于描述方程或系统：

刚性方程（Stiff Equation） 是指数值解只有在时间间隔很小时才会稳定的微分方程。目前很难去精确地去定义哪些微分方程是刚性方程，若此方程式中包含使其快速变动的项，则其为刚性方程。

刚性系统（Stiff System） 通常指由刚性方程描述的系统。刚性系统的动态响应具有多重时间尺度，包括快速衰减或剧烈变化的部分和慢速变化的部分，系统状态变量的变化速率存在显著差异。

刚性系统的特点

刚性系统的特点主要体现在它的动态行为和数值求解的困难性上，尤其是在系统中存在多个不同的时间尺度时表现明显。

刚性系统一般有以下几个特点：

多时间尺度：刚性系统通常有一些状态变量变化得非常快，而另一些变量则变化得很慢。这种差异导致在数值求解时，使用标准方法（如欧拉法或普通的 Runge-Kutta 方法）可能会要求非常小的时间步长来保持数值稳定性，从而增加计算成本。
数值稳定性问题：在刚性系统中，如果使用较大的步长进行仿真，数值解可能会变得不稳定或发散。为了避免这种情况，常规的求解器会自动减小步长以维持稳定性，但这会导致仿真速度显著降低。
需要隐式求解器： 由于刚性系统的数值不稳定性问题，显式求解方法通常效率低或不适用。隐式求解器则更适合刚性系统，因为它们在大步长下仍然能够保持数值稳定性，并且不会因时间尺度差异产生发散问题。
快速衰减的特征值： 刚性系统的线性化方程中往往包含非常大的负特征值（特征值实部非常大且为负），对应于系统的快速衰减成分。这些大负特征值会使显式方法需要使用异常小的时间步长来维持数值稳定。
物理意义上的刚性： 刚性的物理系统有一个或多个在通常意义上表现为“刚性”的组件，例如具有大弹性系数的弹簧、准不可压缩流体和低电感。此类系统通常在其某些组件或模式中表现出高频振荡。
刚性与系统稳定性： 刚性系统通常包含快速衰减的模式，这意味着系统在物理层面上可能是稳定的（快速衰减至稳态），但从数值求解的角度，刚性的使显式方法在求解中表现为不稳定，进而需要隐式方法来解决。

如何判断系统是否为刚性系统

判断一个系统是否是刚性必须要考虑系统的特征值，刚性系统的所有特征值的实部均为负数，并且其中特征值实部绝对值中，最大和最小的比值大于1，该比值也称为刚性比。

刚性比（stiffness ratio）是衡量系统刚性程度的重要指标，定义为：

刚性比最大特征值最小特征值

因此，刚性比越大，系统越刚性，数值求解越困难。

刚性问题的求解

在计算刚性问题时，常见的方法可以分为刚性方法和非刚性方法， 这两种方法的核心区别在于对问题中的刚性约束的处理方式。

刚性方法能够更高效处理刚性问题，同时有着更好的数值稳定性。

非刚性方法（大多是显式方法）

刚性方法（通常是隐式方法）

非刚性方法一般为显式方法（explicit method） ，其原理是利用上一步和当前步的结果计算下一步的计算结果，用公式可以表示为：

显式方法的特点是：

基于过去的信息来计算下一步
计算量较小
稳定性较差
非常适合于非刚性系统

刚性方法一般为隐式方法（inplicit method）， 隐式方法在计算下一步时，需要同时满足当前和未来的状态，即通过迭代方法来求解下一步，使其在给定的时间步内达到一致性。用公式可以表示为：

在上式中，出现在方程的右侧，因此需要通过迭代求解出下一步。

隐式方法的特点是：

能自洽地求解下一步
计算量较大
稳定性强
更适合刚性系统，特别是对于物理系统模型的求解。

需要注意的是：

刚性问题是一个效率问题，如果不关心计算效率，就不需要关心刚性。

非刚性方法也可以解决刚性问题，只是需要更长的时间。

刚性的影响

前文提到，系统或方程的刚性主要会带来数值稳定性与计算效率的问题

数值稳定性

这里以一个简单的案例进行解释：

这里有一个微分方程组：

可以看到两个方程的时间尺度相差较大。

其特征值为：

前文提到刚性比（stiffness ratio）是衡量系统刚性的重要指标，根据其定义可得该系统的刚性比为：

刚性比 = |最大特征值| / |最小特征值| = 1000/1 = 1000

因此该系统属于刚性系统。

其解析解为：

其中：

是一个指数衰减函数，从初始值开始逐渐衰减。

由两部分组成，其中在初始时刻占主导地位，但其是一个快速衰减的项，当稍微增大时（比如），就迅速趋近于 0。

所以，和的解只会在初始阶段(0~0.01s)有所不同，但随着时间的推移，它们的值会逐渐趋于相同。

在初始阶段会在极短的时间内衰减至1，如此迅速的变化会产生极大的刚性，给数值求解的稳定性带来了困难。

如果使用显式方法进行计算，对于刚性问题，显式欧拉方法的数值稳定性条件是由特征值的模所决定的。

一般来说，步长需要满足如下的稳定性条件：

其中是特征值中模最大的值。

在这个问题中，特征值为和，模最大值为 1000。因此步长的最大值为：

因此，使用显式欧拉方法时，步长需要小于 0.002s才能保证数值稳定性。

我们在Simulink中搭建该系统模型进行验证：

选择非刚性方法的显式欧拉求解器ode1。

当步长取0.001s时（小于最大步长），计算结果为：

显式欧拉法步长0.001s

案例的解析解为：

解析解

通过对比解析解与数值解可以看到非刚性显式方法在使用正确的步长时可以正确地求解，特别是y2由初始值2快速衰减至1的过程。

但当步长为0.002s，计算结果为：

显式欧拉法步长0.002s

可以看到y2由于其方程的高刚性，导致了数值解在[-1,2]的范围内高频振荡。

通过观察[0,0.02s]的图像，可以看到振荡的细节，其数值解以0.002s的步长反复振荡。

步长0.001s, [0,0.02s]

造成震荡的原因是：

首先，初始值：y₁(0) = 1, y₂(0) = 2，步长：h = 0.002

根据显式欧拉法的公式计算y₂：

第一步：

第二步：

第三步：

根据显式欧拉法计算得到前三步的数值解：

数值振荡的振幅由以下因素决定：

初始值
步长（h）
方程中的系数（999和-1000）

基于上述结论，可以证明随着的衰减，本案例中的震荡范围为[-1，2]。

如果使用更高的步长，例如h = 0.01s，数值不稳定的现象会更加明显，求解将很难进行下去。

y2数值不稳定

求解异常

计算效率

前文提到，显式方法在求解刚性问题时稳定性较差，需要较小的步长。隐式方法稳定性较高，可以使用较大的步长，但单个步长的计算量会高于显式方法。

因此，虽然显式方法的总步数更多，但单个步长的计算量较低，而隐式方法总步数更少，尽管单步的计算成本较高，但在大部分情况下整体效率更高。

这里使用隐式方法求解前文的案例

选择隐式后向欧拉方法ode1be，后向欧拉法的公式为：

后向欧拉法可以选择较大的步长，即使h = 1s时，虽然牺牲了一定的精度，但求解依旧保持了稳定。

隐式欧拉法步长1s

当步长选择h = 0.02s时，对比解析解其精度已经可以满足需求。

隐式欧拉法步长0.02s

相较于步长h = 0.001s才能稳定求解的显式欧拉法，隐式欧拉法使用h = 0.02s的步长就能达到相当的精度。

借助Simulink的Simulink Profiler可知在求解阶段simulationPhase，隐式方法比显式方法有更高的计算效率。

显式方法

隐式方法

借助Simulink的Solver Profiler求解器探查工具，可知显式方法步长更多，隐式方法步长少但需要更新雅可比矩阵，但整体计算效率更高（运行时间）。

显式方法

隐式方法

实际情况下，在面对刚性问题选择求解方法时需要根据问题的实际情况与计算资源进行权衡，需要考虑的因素有：

问题的类型： 如果问题是高度动态的，或者需要捕捉细节丰富的瞬态现象，则选择显式求解器，更多的步数能获得较好的时间分辨率。如果问题是趋于稳态或者涉及到长时间尺度的行为，选择隐式求解器更高效，因为它在大时间步长下可以更稳定地求解，减少计算步数。
所需精度： 显式方法的高精度依赖于使用小时间步长来捕捉快速变化的现象，因此在一些精度要求高的瞬态问题中，它的优势较为明显。隐式方法可以通过调整步长在精度和效率之间找到平衡。
计算资源： 如果有大量并行计算资源，显式方法通常更容易并行化，因为每一步的计算独立性较强，可以有效地利用多核或多节点的计算资源。
问题的非线性程度： 对于强非线性问题，显式方法不需要迭代求解每一步的线性化方程组，往往更容易处理这些非线性特征，实现和调试更简单。隐式方法在处理非线性问题时，也可以通过非线性迭代方法来稳定求解，但实现起来更复杂。
代码复杂性： 如果简单性是首要考虑因素，显式方法更适合追求简单性和快速原型的场合。
内存限制： 如果内存或存储空间有限制，显式方法可能更有优势，隐式方法可能需要大量的存储空间用于存储系数矩阵和中间变量