刚性检测

概述"NonstiffTest" 方法选项
引言示例
线性稳定性选项总结
"StiffnessTest" 方法选项

概述

许多微分方程表现出某些形式的刚性，从而限制步长大小，进而限制了显式求解方法的有效性.

多年来，大量隐式方法得到开发，以克服这个问题.

对于同样的步长，由于与内在的线性代数相关的开销，隐式方法的效率可能会显著低于显式方法.

由于在某些区域，隐式方法可以使用显著增加的步长，该成本可以被抵消.

我们作了多次尝试来提供用户友好的代码，以便自动尝试在运行时检测刚性，并在必要时在适当的方法之间转换.    

这里对一些已经提出的自动为代码配备刚性检测设备的策略作一概述.

特别关注的是一个矩阵的主导特征值的估计问题，以描述刚性检测如何在 NDSolve 里实现.

数值实例展现了该策略的有效性.

初始化

加载具有预定义示例和功用函数的程序包：

引言

考虑初值问题：

的数值解.刚性是问题、解法、初始条件和局部误差容差的一个结合体.

由于对可以采取的步长的限制，刚性限制了显式解法的有效性.

刚性出现在许多实际系统中，以及在由直线法求解的偏微分方程的数值解中.

示例

范德波尔振荡器是一个非保守的非线性阻尼振荡器，它是一个常微分方程：

的刚性系统的例子. 其中 ；

考虑初始条件：

并且在区间 上求解.

方法 "StiffnessSwitching" 默认时使用一对插值方法：

解

这里从一个程序包加载该问题：
使用非刚性方法数值求解该系统：
使用一种当刚性发生时进行转换的方法求解该系统：
下面是两个解组成部分的图形. 尖峰（蓝色）在幅度上扩大至大约450，并且已经被裁剪：

刚性往往发生在快速瞬变后的区域.

这里对采取的步长和时间的关系画出图形：

问题是，当解快速变化时，几乎没有用刚性求解器的必要，因为局部准确性是主要问题.

为了提高效率，如果该方法能够自动检测出那些局部准备性（而不是稳定性）更为重要的区域，那将是很有用的.

线性稳定性

线性稳定性理论源于 Dahlquist 标量线性测试方程：

的研究，作为研究初值问题 (1) 的一个简化模型.

稳定性是通过分析应用到 (2) 以得到下式的方法来描述的：

这里 ， 是（有理）稳定性函数.

绝对稳定性的边界由考虑如下区域获得：

显式欧拉法

显式或者向前欧拉法：

用于 (3) 给出：

阴影区域代表不稳定性，其中 >1.

线性稳定性边界往往取为与负实轴的交点.

对于显式欧拉方法，线性稳定性边界 .

对于一个特征值 ，线性稳定性的要求表示步长需要满足 ，这是一个非常温和的限制.

然而，对于一个特征值 ，线性稳定性的要求表示步长需要满足 ，这是一个非常严格的限制.

示例

这个例子表明当使用显式 Runge-Kutta 方法求解一个刚性系统时，刚性对步长序列的影响.

该系统对于一个化学反应建模：
该系统通过禁用内置的刚性检测求解：
当达到稳定边界时，步长序列开始振荡：

隐式欧拉法

隐式欧拉法或者后退欧拉法：

用于 (4) 给出：

该方法对于整个左侧半平面是无条件稳定的.

这意味着若要保持稳定性，不再需要有对步长的限制.

缺点是在每个积分步骤都要求解一个隐式方程组.

类型不敏感（Type Insensitivity）

一个对类型不敏感的求解器在每个步骤对刚性进行确认并有效地作出反应，因此对于问题的类型（有可能是变化的）是不敏感的.

这种类型的求解器中最完善的之一是 LSODA [H83], [P83].

LSODA 的后代求解器如 CVODE 不再纳入刚性检测设备. 原因是 LSODA 使用范数界限来估计主导特征值，而这些界限可以是相当不准确的，这一点将在后面的讨论中看到.

低阶的A(α)-稳定的 BDF 方法意味着 LSODA 和 CVODE 不适合求解具有高准确度的，或者主导特征值具有大的虚部的系统. 其他的方法，如那些基于线性隐式的方案的外推法就不会有这些问题.    

关于刚性检测的大部分工作是在80年代和90年代使用独立的 FORTRAN 代码实现的.

此后，出现了新的线性代数技术和有效的软件，这些都可在 Wolfram 语言中获得.

刚性可以是一个短暂的现象，因此检测非刚性是同样重要的 [S77], [B90].

"StiffnessTest" 方法选项

有好几种方法可用于从非刚性求解器到刚性求解器的转变.

直接估计

检测刚性的一个方便的方法就是直接估计该问题的雅克比 的主导特征值（见 [S77], [P83], [S83], [S84a], [S84c], [R87] 和 [HW96]）.    

这一个估计往往可以作为数值积分的副产品加以利用，所以代价合理地低.

如果 代表对应于雅可比的主导特征值的特征向量的一个近似，同时 足够小，那么根据中值定理，主导特征值的一个很好的近似是：

Richardson 外推法提供了产生这种形式的量的一个优化序列，如某些显式 Runge–Kutta方法一样

成本最多是两个函数计算，但是往往其中至少有一个是数值积分的副产品，因此它的代价合理地低.

用 表示线性稳定性边界——它是线性稳定性区域和负实轴的交集.

积 给出一个可以和一个方法的线性稳定性边界相比较的估计，以便检测刚性：

其中 是一个安全因子.

描述

方法 "DoubleStep"、"Extrapolation" 和 "ExplicitRungeKutta" 具有选项 "StiffnessTest"，它可以用来识别该方法在应用指定的 AccuracyGoal 和 PrecisionGoal 容差时对于一个给定的问题是否是刚性的.

方法选项 "StiffnessTest" 本身接受许多选项来执行 (5) 的一个弱形，允许测试失败指定数目的次数.

这样做的原因是有些问题只在某些区域有轻度的刚性，而一个显式积分方法可能仍然是有效的.

"NonstiffTest" 方法选项

"StiffnessSwitching" 方法具有选项 "NonstiffTest"，它用来把一个刚性方法转变回一个非刚性方法.

对于选项 "NonstiffTest" 允许下列设置：

转化到一个非刚性求解器

这里采用一个独立于刚性方法的方法.

给定雅可比 J（或者一个近似）计算如下之一：

范数界限（Norm Bound）：

谱半径：ρ(J)max

主导特征值 .

许多线性代数技术专注于高精地解决一个单一的问题.

对于刚性检测，具有精确到一或二位数字的解的一系列问题就足够了.

对于一个数值离散化

考虑在某（一些）子区间上的一系列 个矩阵

一系列相连的矩阵的谱经常慢慢地一步一步改变.

我们的目标是找到一种估算 矩阵系列 的主导特征值（的边界），以满足：

NormBound

获取主导特征值的界限的一个简单有效的技术是使用雅可比的范数 ， 通常 或者 .

该方法具有复杂度 ，它少于在刚性求解器中采取的工作.

下面是 LSODA 所使用的方法.

示例

下面雅可比矩阵出现在使用刚性求解器的 Van der Pol 系统的数值解中：
基于范数的界限对谱半径的高估超过一个数量级：

直接特征值计算

对于小型问题 () 直接计算主导特征值可能更有效.

幂法

Shampine 提出了采用幂法来估计雅可比的主导特征值 [S91].

幂法也许不是一个很备受尊重的方法，但是由于它在 Google 页面排名中的使用，很多人已经重新对其表示兴趣.

幂法可以用于当

描述

给定一个初始向量 ，计算

用 Rayleigh 商来计算主导特征值的近似：

在实践中，近似特征向量在每一步被归一化：

属性

幂法线性收敛，速率为：

这可能会很慢.

特别是，当应用于具有复共轭的主导特征值对时，该方法不收敛.

推广

幂法可以通过适应性修正来克服等模（equimodular）特征值问题（如 NAPACK）.

但是，这样的修正一般不处理集群特征值的收敛速度慢的问题.

有两种主要的方法来推广幂法:

子空间迭代

子空间（或同时）迭代通过每一步作用于 m 个向量推广了幂法中的思想.

从一个正交向量集 开始，通常 :

形成与矩阵 的积：

为了防止所有向量收敛到 的同一主导特征向量 的倍数，它们被正交归一化：

比起矩阵乘积来说，正交归一化的步骤是比较昂贵的.

Rayleigh–Ritz 投影

输入：矩阵 和向量 的一个标准正交集合

收敛性

子空间（或者同时）迭代通过在每一步作用于 个向量，推广了幂法中的思想.

SRRIT 线性收敛，速度为：

特别地，主导特征值的收敛速度为：

因此，采用比如 或者更多是有益的，即使我们只对主导特征值感兴趣.

误差控制

在连续近似主导特征值上的相对误差的测试为：

这是不够的，因为当收敛速度很慢时，它仍可以被满足.

如果 或者 那么 的第 列不是唯一确定的.

用于 SRRIT 的残差测试为：

其中 ， 是 的第 列，而 是 的第 列.

这是有利的，因为它适用于等模（equimodular）特征值.

因为有序 Schur 分解的使用，要对上三角矩阵 的第一列的位置进行测试.

实施

子空间迭代的实现方法有好几种.

KrylovIteration

给定一个 矩阵 ，其列 组成了一个给定子空间 的正交基：

Rayleigh Ritz 程序包含计算 以及求解相关的特征值问题 .

原问题的近似特征对 , 满足 和 ，它们被称为 Ritz 值和 Ritz 向量.

当子空间 接近 的不变子空间时，该过程效果最佳.

当 等于与矩阵 和给定初始向量 相关的 Krylov 子空间时，这个过程有效：

描述

Arnoldi 方法是一种基于 Krylov 的投影算法，它计算 Krylov 子空间的一个正交基，并且产生一个投影的 矩阵 其中 .

输入：矩阵 ，步骤数目 ，模为1的初始向量 .

输出： 其中

在 Arnoldi 的情况下， 具有一个未约化（unreduced）的上黑森贝格形式（ upper Hessenberg form），即具有一个额外的非零次对角的上三角.

正交化通常是通过格拉姆-施密特程序的方法来完成的.

由该算法计算的量满足：

残差 给出了一个对不变子空间的近似的指示，相关的范式 表明了计算 Ritz 对的精度：

重启

如果初始向量 在预期特征值的方向是丰富的话，Ritz 对的收敛速度很快的.

当情况并非如此时，则需要一个重新启动的策略，以避免工作量和内存的过快增长.

重启策略有好几种，特别地有：

实现

有许多软件实现可用，尤其是：

Automatic 策略

"Automatic" 设置使用这些方法的一个组合，情形如下.

步骤拒绝

对计算时间的缓存确保了主导特征值估计不在被拒绝的步骤重新计算.

对被拒绝的步骤也要进行刚性检测，这是由于：

迭代方法选项

"NonstiffTest" 的迭代方法具有可被修改的选项：

默认容差的目的是只有一个正确数位，但往往得到更准确的值——特别是在连续的步骤中几个成功的迭代之后.

限制迭代数目的默认值是：

延迟和切换

有一点很重要，就是将某种形式的延迟结合进来，以避免一种 "StiffnessSwitching" 方法不断尝试在刚性和非刚性方法之间切换的循环.

为此目的，我们使用 "StiffnessTest" 和 "NonstiffTest" 的选项 "MaxRepetitions" 和 "SafetyFactor".

默认设置允许切换是相当被动的，这对于单步积分方法是恰当的.

示例

Van der Pol

选择一个实例系统：

StiffnessTest

该系统用一个给定的方法和 "StiffnessTest" 的默认选项设置被成功地积分：
当刚性测试条件不满足时，一个较长的积分将中止，并发出提示信息：
使用单位安全系数，并指定只允许有一个刚性失败，这样做实际上给出了一个严格的测试. 该项指定使用了嵌套的方法选择句法：

NonstiffTest

对于这样的一个小系统，采用直接特征值计算.

这个例子是一个很好的测试，体现了整体刚性切换架构的行为正如预期.

建立一个函数来监控刚性方法和非刚性方法之间的转换以及所采用的步长. 刚性方法和非刚性方法的数据通过使用 Sow 的一个不同标签放置在了不同的列表中：
求解系统，并且收集用于方法转换的数据：
画出使用显式求解器（蓝色）和隐式求解器（红色）所采用的步长的图形：
计算采用的非刚性和刚性步的数目（包括拒绝步）：

CUSP

关于神经冲动机制的 cusp catastrophe 模型 [Z72]：

与 Van der Pol 振荡器相结合产生了 CUSP 系统 [HW96]：

其中

而 并且 .

利用直线法对扩散项进行离散化从而得到一个维度为 的常微分方程组.     

与 Van der Pol 系统不同，因为问题的大小，我们采用迭代方法来进行特征值估计.

步长和阶数选择

选择要求解的问题：
建立一个函数来监控使用的方法类型和步长. 另外，方法的阶数作为一个工具提示条被包含在内：
收集当积分进行时方法阶数（order）的数据：
画出使用显式求解器（蓝色）和隐式求解器（红色）所采用的步长. 一个工具提示条显示了在每步的方法阶数：
计算所采用的非刚性和刚性步的总数（包括拒绝步）：

雅可比实例

定义一个函数收集前几个雅可比（Jacobian）矩阵：

一个到刚性方法的转换出现在 0.00113425 附近，而非刚性的第一个测试出现在下一步 .

雅可比 的图形展示：
定义一个函数来计算并且显示 , , … 等等的前几个特征值，以及范数界限：

范数界限在这个例子中非常明显.

Korteweg–deVries

Korteweg–deVries 偏微分方程是在浅水表面的波浪的数学模型：

我们考虑如下边界条件：

并且在区间 上求解.

采用直线法的离散化用来形成由 192 个常微分方程组成的方程组.

步长

选择要求解的问题：
在 LSODA 中使用的向后差分公式方法在求解这一问题时遇到了困难：
一个图形显示了步长快速减小：
相反，StiffnessSwitching 立即切换到使用线性隐式欧拉方法，它需要很少的积分步骤：

一旦在积分开始的时候选择了刚性方法，插值法从来不会切换回非刚性求解器.

因此，这是一个非刚性检测的最差情况例子的形式.

尽管如此，使用子空间迭代的成本仅是总积分时间的百分之几.

在关闭到非刚性方法的转换时，计算所采用的时间：

雅可比实例

采用前面定义的监控函数，收集前几个雅可比矩阵的数据：
初始雅可比 的图形显示：
计算并且显示 , , … 的前几个特征值，以及范数界限：

范数界限略被高估，但更重要的是，他们没有提供实部和虚部的相对大小的显示：

选项总结

StiffnessTest

选项名
默认值
"MaxRepetitions"{3,5}指定刚性测试 (6) 允许的连续失败的最大次数和失败的总次数的最大值
"SafetyFactor"指定在刚性测试 (7) 的右侧所使用的安全系数

方法选项 "StiffnessTest" 的选项.

NonstiffTest

选项名
默认值
"MaxRepetitions"{2,∞}指定刚性测试 (8) 允许的连续失败的最大次数和失败的总次数的最大值
"SafetyFactor"指定在刚性测试(9)的右侧所使用的安全系数

方法选项 "NonstiffTest" 的选项.