Wolfram Mathematica Documentation Center
DOCUMENTATION CENTER SEARCH
New to Mathematica? Find your learning path »
虚拟全书 > 数学和算法 > Mathematica 中的高等微分方程的数值求解 > 刚性检测 >
MATHEMATICA 教程

刚性检测

概述

许多微分方程表现出某些形式的刚性，从而限制步长大小，进而限制了显式求解方法的有效性.
多年来，大量隐式方法得到开发，以克服这个问题.
对于同样的步长，由于与内在的线性代数相关的开销，隐式方法的效率可能会显著低于显式方法.
由于在某些区域，隐式方法可以使用显著增加的步长，该成本可以被抵消.
我们作了多次尝试来提供用户友好的代码，以便自动尝试在运行时检测刚性，并在必要时在适当的方法之间转换.    
这里对一些已经提出的自动为代码配备刚性检测设备的策略作一概述.
特别关注的是一个矩阵的主导特征值的估计问题，以描述刚性检测如何在 NDSolve 里实现.
数值实例展现了该策略的有效性.

初始化

加载具有预定义示例和功用函数的程序包.
In[59]:=
Click for copyable input

引言

考虑初值问题：
的数值解.
刚性是问题、解法、初始条件和局部误差容差的一个结合体.
由于对可以采取的步长的限制，刚性限制了显式解法的有效性.
刚性出现在许多实际系统中，以及在由直线法求解的偏微分方程的数值解中.

示例

范德波尔振荡器是一个非保守的非线性阻尼振荡器，它是一个常微分方程：
的刚性系统的例子. 其中 ；
考虑初始条件：
并且在区间 上求解.
方法 默认时使用一对插值方法：
  • 显式修改中点（Explicit modified midpoint，即 Gragg 平滑），双谐波序列 2，4，6 ...
  • 线性隐式欧拉，次谐波序列 2，3，4 ...

解

这里从一个程序包加载该问题.
In[10]:=
Click for copyable input
使用非刚性方法数值求解该系统.
In[11]:=
Click for copyable input
使用一种当刚性发生时进行转换的方法求解该系统.
In[12]:=
Click for copyable input
下面是两个解组成部分的图形. 尖峰（蓝色）在幅度上扩大至大约450，并且已经被裁剪.
In[13]:=
Click for copyable input
Out[13]=
刚性往往发生在快速瞬变后的区域.
这里对采取的步长和时间的关系画出图形.
In[14]:=
Click for copyable input
Out[14]=
问题是，当解快速变化时，几乎没有用刚性求解器的必要，因为局部准确性是主要问题.
为了提高效率，如果该方法能够自动检测出那些局部准备性（而不是稳定性）更为重要的区域，那将是很有用的.

线性稳定性

线性稳定性理论源于 Dahlquist 标量线性测试方程：
的研究，作为研究初值问题 (1) 的一个简化模型.
稳定性是通过分析应用到 (2) 以得到下式的方法来描述的：
这里 ， 是（有理）稳定性函数.
绝对稳定性的边界由考虑如下区域获得：

显式欧拉法

显式或者向前欧拉法：
用于 (3) 给出：
阴影区域代表不稳定性，其中 .
In[15]:=
Click for copyable input
Out[15]=
线性稳定性边界往往取为与负实轴的交点.
对于显式欧拉方法，线性稳定性边界 .
对于一个特征值 ，线性稳定性的要求表示步长需要满足 ，这是一个非常温和的限制.
然而，对于一个特征值 ，线性稳定性的要求表示步长需要满足 ，这是一个非常严格的限制.

示例

这个例子表明当使用显式 Runge-Kutta 方法求解一个刚性系统时，刚性对步长序列的影响.
该系统对于一个化学反应建模.
In[16]:=
Click for copyable input
该系统通过禁用内置的刚性检测求解.
In[17]:=
Click for copyable input
当达到稳定边界时，步长序列开始振荡.
In[18]:=
Click for copyable input
Out[18]=
  • 大量的步拒绝（step rejections）常常会对性能产生负面影响.
    • 大量步骤的使用对计算结果的准确性产生不利影响.
    内置检测对刚性出现时间的定位表现出色.
    In[19]:=
    Click for copyable input

隐式欧拉法

隐式欧拉法或者后退欧拉法：
用于 (4) 给出：
该方法对于整个左侧半平面是无条件稳定的.
In[20]:=
Click for copyable input
Out[20]=
这意味着若要保持稳定性，不再需要有对步长的限制.
缺点是在每个积分步骤都要求解一个隐式方程组.

类型不敏感（Type Insensitivity）

一个对类型不敏感的求解器在每个步骤对刚性进行确认并有效地作出反应，因此对于问题的类型（有可能是变化的）是不敏感的.
这种类型的求解器中最完善的之一是 LSODA [H83], [P83].
LSODA 的后代求解器如 CVODE 不再纳入刚性检测设备. 原因是 LSODA 使用范数界限来估计主导特征值，而这些界限可以是相当不准确的，这一点将在后面的讨论中看到.
低阶的A()-稳定的 BDF 方法意味着 LSODA 和 CVODE 不适合求解具有高准确度的，或者主导特征值具有大的虚部的系统. 其他的方法，如那些基于线性隐式的方案的外推法就不会有这些问题.    
关于刚性检测的大部分工作是在80年代和90年代使用独立的 FORTRAN 代码实现的.
此后，出现了新的线性代数技术和有效的软件，这些都可在 Mathematica 中获得.
刚性可以是一个短暂的现象，因此检测非刚性是同样重要的 [S77], [B90].

"StiffnessTest" 方法选项

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

直接估计

检测刚性的一个方便的方法就是直接估计该问题的雅克比 的主导特征值（见 [S77], [P83], [S83], [S84a], [S84c], [R87] 和 [HW96]）.    
这一个估计往往可以作为数值积分的副产品加以利用，所以代价合理地低.
如果 代表对应于雅可比的主导特征值的特征向量的一个近似，同时 足够小，那么根据中值定理，主导特征值的一个很好的近似是：
Richardson 外推法提供了产生这种形式的量的一个优化序列，如某些显式 Runge-Kutta方法一样
成本最多是两个函数计算，但是往往其中至少有一个是数值积分的副产品，因此它的代价合理地低.
用 表示线性稳定性边界——它是线性稳定性区域和负实轴的交集.
积 给出一个可以和一个方法的线性稳定性边界相比较的估计，以便检测刚性：
其中 是一个安全因子.

描述

方法 、 和 具有选项 ，它可以用来识别该方法在应用指定的 AccuracyGoal 和 PrecisionGoal 容差时对于一个给定的问题是否是刚性的.
方法选项 本身接受许多选项来执行 (5) 的一个弱形，允许测试失败指定数目的次数.
这样做的原因是有些问题只在某些区域有轻度的刚性，而一个显式积分方法可能仍然是有效的.

"NonstiffTest" 方法选项

方法具有选项 ，它用来把一个刚性方法转变回一个非刚性方法.
对于选项 允许下列设置：

转化到一个非刚性求解器

这里采用一个独立于刚性方法的方法.
给定雅可比 J（或者一个近似）计算如下之一：
范数界限（Norm Bound）：
谱半径：
主导特征值 .
许多线性代数技术专注于高精地解决一个单一的问题.
对于刚性检测，具有精确到一或二位数字的解的一系列问题就足够了.
对于一个数值离散化
考虑在某（一些）子区间上的一系列 个矩阵
一系列相连的矩阵的谱经常慢慢地一步一步改变.
我们的目标是找到一种估算矩阵系列的主导特征值（的边界），以满足：
  • 在刚性求解器中的每一步，成本低于在线性代数中实施的工作.
  • 将求解器按步求解的性质考虑在内.

NormBound

获取主导特征值的界限的一个简单有效的技术是使用雅可比的范数 ， 通常 或者 .
该方法具有复杂度 ，它少于在刚性求解器中采取的工作.
下面是 LSODA 所使用的方法.
  • 稠密矩阵的范数界限高估了，而在维度增加时，界限变得更不理想.
    • 范数界限对于具有相当大维度的稀疏矩阵或者带状矩阵可能是很紧张的.
    选项 的设置 计算 和 并且返回这两个值中较小的一个.

示例

下面雅可比矩阵出现在使用刚性求解器的 van der Pol 系统的数值解中.
In[21]:=
Click for copyable input
基于范数的界限对谱半径的高估超过一个数量级.
In[22]:=
Click for copyable input
Out[22]=

直接特征值计算

对于小型问题 () 直接计算主导特征值可能更有效.
  • 厄米 （Hermitian） 矩阵采用 LAPACK 函数 xgeev
    • 一般矩阵使用 LAPACK 函数 xsyevr
    选项 的设置 采用与 Eigenvalues 相同的 LAPACK 程序计算 的主导特征值.
    对于较大的问题，直接计算特征值的代价是 ，当与在刚性求解器中进行的线性代数的代价相比时，变得不可接受.
    为此，我们已经实现了许多迭代方案. 这些实际是通过在一个较小的子空间中对主导特征空间进行近似，以及对于较小的问题使用稠密特征值方法来工作的.

幂法

Shampine 提出了采用幂法来估计雅可比的主导特征值 [S91].
幂法也许不是一个很备受尊重的方法，但是由于它在 Google 页面排名中的使用，很多人已经重新对其表示兴趣.
幂法可以用于当
  • 具有 n 个线性独立的特征向量（可对角化的）.
  • 特征值可以按大小排序，如 .
  • 是 的主导特征值.

描述

给定一个初始向量 ，计算
用 Rayleigh 商来计算主导特征值的近似：
在实践中，近似特征向量在每一步被归一化：

属性

幂法线性收敛，速率为：
这可能会很慢.
特别是，当应用于具有复共轭的主导特征值对时，该方法不收敛.

推广

幂法可以通过适应性修正来克服等模（equimodular）特征值问题（如 NAPACK）.
但是，这样的修正一般不处理集群特征值的收敛速度慢的问题.
有两种主要的方法来推广幂法:
  • 对于小到中维度的子空间迭代
    • 对于大维度的 Arnoldi 迭代
    虽然方法的工作方式完全不同，但是仍然有一些核心部分可以被共享和优化.
    子空间和 Krylov 迭代的成本为 步操作.    
    它们把一个 矩阵投影到一个 矩阵，其中 .
    小矩阵表示主导特征空间，近似中采用稠密特征值程序.

子空间迭代

子空间（或同时）迭代通过每一步作用于 m 个向量推广了幂法中的思想.
从一个正交向量集 开始，通常 :
形成与矩阵 的积：
为了防止所有向量收敛到 的同一主导特征向量 的倍数，它们被正交归一化：
比起矩阵乘积来说，正交归一化的步骤是比较昂贵的.

Rayleigh-Ritz 投影

输入：矩阵 A 和向量 V 的一个标准正交集合
  • 计算 Rayleigh 商
    • 计算 Schur 分解
    矩阵 S 具有小的维度 .
    请注意，当 时，使用一个准上三角矩阵 ，Schur 分解可以实数运算来计算.    

收敛性

子空间（或者同时）迭代通过在每一步作用于 个向量，推广了幂法中的思想.
SRRIT 线性收敛，速度为：
特别地，主导特征值的收敛速度为：
因此，采用比如 或者更多是有益的，即使我们只对主导特征值感兴趣.

误差控制

在连续近似主导特征值上的相对误差的测试为：
这是不够的，因为当收敛速度很慢时，它仍可以被满足.
如果 或者 那么 的第 列不是唯一确定的.
用于 SRRIT 的残差测试为：
其中 ， 是 的第 列，而 是 的第 列.
这是有利的，因为它适用于等模（equimodular）特征值.
因为有序 Schur 分解的使用，要对上三角矩阵 的第一列的位置进行测试.

实施

子空间迭代的实现方法有好几种.
  • 采用切比雪夫（Chebyshev） 加速的子空间迭代 [S84b], [DS93]
    在 中所用的实现方法是以：
    • Schur Rayleigh-Ritz 迭代 [BS97]
    为基础的.
    "SRRIT 的一个吸引人的特点是它显示了单调一致性，即当收敛容差减小时，计算得的残差的大小也减小" [LS96].    
    SRRIT 用到了有序 Schur 分解，其中具有最大模的特征值出现在左上角的项中.
    的形成用到了重新正交归一化的修正的 Gram-Schmidt，这比Householder变换速度更快.
    在积分时间 t_i 的近似主导子空间 用于启动在下一个积分步骤 t_(i+1) 时的迭代：

KrylovIteration

给定一个 矩阵 ，其列 组成了一个给定子空间 的正交基：
Rayleigh Ritz 程序包含计算 以及求解相关的特征值问题 .
原问题的近似特征对 , 满足 和 ，它们被称为 Ritz 值和 Ritz 向量.
当子空间 接近 的不变子空间时，该过程效果最佳.
当 等于与矩阵 和给定初始向量 相关的 Krylov 子空间时，这个过程有效：

描述

Arnoldi 方法是一种基于 Krylov 的投影算法，它计算 Krylov 子空间的一个正交基，并且产生一个投影的 矩阵 其中 .
输入：矩阵 ，步骤数目 ，模为1的初始向量 .
输出： 其中
在 Arnoldi 的情况下， 具有一个未约化（unreduced）的上黑森贝格形式（ upper Hessenberg form），即具有一个额外的非零次对角的上三角.
正交化通常是通过格拉姆-施密特程序的方法来完成的.
由该算法计算的量满足：
残差 给出了一个对不变子空间的近似的指示，相关的范式 表明了计算 Ritz 对的精度：

重启

如果初始向量 在预期特征值的方向是丰富的话，Ritz 对的收敛速度很快的.
当情况并非如此时，则需要一个重新启动的策略，以避免工作量和内存的过快增长.
重启策略有好几种，特别地有：
  • 显式重启 —— 一个新的启动向量是 Ritz 向量的子集的线性组合.
    • 隐式重启 —— 从与隐式位移 QR 算法相结合的 Arnoldi 过程中形成的一个新的启动向量.
    显式重启的实现相对简单，但是隐式重启更有效，因为它保留了更大问题的相关特征信息. 然而，隐式重启以数值稳定的方式实现起来会很困难.
    另一种方法更容易实现，但是可以达到与隐式重启相同的效果，它是 Krylov-Schur方法 [S01].

实现

有许多软件实现可用，尤其是：
    在 中的实现是基于Krylov-Schur 迭代的.

Automatic 策略

设置使用这些方法的一个组合，情形如下.
  • 对于 采用直接特征值计算. 或者使用 或者 ，取决于哪种方法是主动（active）的.
  • 对于 采用子空间迭代，默认的基的大小为 . 如果该方法成功，那么得到的基用来在下个积分步骤启动该方法.
  • 如果子空间迭代在 max_(si) 个迭代后没有收敛，那么主导向量将被用来启动默认基大小为 的 Krylov 方法. 接下来的积分步骤使用 Krylov 方法，从上一步得到的向量开始.
  • 如果 Krylov 迭代在 个迭代后没有收敛，那么范数界限将被用于当前步骤. 下一个积分步骤将继续尝试 Krylov 迭代.
  • 由于它们的代价如此低，当使用子空间或者 Krylov 迭代时，以及绝对值中较小的被使用时，范数界限总是被计算.

步骤拒绝

对计算时间的缓存确保了主导特征值估计不在被拒绝的步骤重新计算.
对被拒绝的步骤也要进行刚性检测，这是由于：
  • 当在稳定性边界工作时，对于非刚性求解器常常出现步骤拒绝.
  • 当求解快速瞬态时，对于刚性求解器常常出现步拒绝.

迭代方法选项

的迭代方法具有可被修改的选项：
In[23]:=
Click for copyable input
Out[23]=
In[24]:=
Click for copyable input
Out[24]=
默认容差的目的是只有一个正确数位，但往往得到更准确的值——特别是在连续的步骤中几个成功的迭代之后.
限制迭代数目的默认值是：
  • 对于子空间迭代 .
    • 对于 Krylov 迭代 .
    如果这些值被设置得太大，那么收敛失败的代价就变得很大.
    在困难的问题中，最好是在步骤之间分享收敛的工作. 由于这些方法实际上改善了上一步的基向量，后续步骤有很大的收敛机会.

延迟和切换

有一点很重要，就是将某种形式的延迟结合进来，以避免一种 方法不断尝试在刚性和非刚性方法之间切换的循环.
为此目的，我们使用 和 的选项 和 .
默认设置允许切换是相当被动的，这对于单步积分方法是恰当的.
  • 对于非刚性方法， 是在每个步的结尾执行的. 当选项 的任何值达到时，出现一个步骤拒绝，然后用一个刚性方法重新计算该步骤.
  • 是抢先执行的. 它是在一个刚性方法采取一个步骤之前， 用上一步的雅可比矩阵来实施的.

示例

Van der Pol

选择一个实例系统.
In[25]:=
Click for copyable input

StiffnessTest

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

NonstiffTest

对于这样的一个小系统，采用直接特征值计算.
这个例子是一个很好的测试，体现了整体刚性切换架构的行为正如预期.
建立一个函数来监控刚性方法和非刚性方法之间的转换以及所采用的步长. 刚性方法和非刚性方法的数据通过使用 Sow 的一个不同标签放置在了不同的列表中.
In[29]:=
Click for copyable input
求解系统，并且收集用于方法转换的数据.
In[31]:=
Click for copyable input
画出使用显式求解器（蓝色）和隐式求解器（红色）所采用的步长的图形.
In[33]:=
Click for copyable input
Out[33]=
计算采用的非刚性和刚性步的数目（包括拒绝步）.
In[34]:=
Click for copyable input
Out[34]=

CUSP

关于神经冲动机制的 cusp catastrophe 模型 [Z72]：
与 van der Pol 振荡器相结合产生了 CUSP 系统 [HW96]：
其中
而 并且 .
利用直线法对扩散项进行离散化从而得到一个维度为 的常微分方程组.     
与 van der Pol 系统不同，因为问题的大小，我们采用迭代方法来进行特征值估计.

步长和阶数选择

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

雅可比实例

定义一个函数来收集前面几个雅可比矩阵.
In[42]:=
Click for copyable input
In[44]:=
Click for copyable input
一个到刚性方法的转换出现在 0.00113425 附近，而非刚性的第一个测试出现在下一步 .
雅可比 的图形展示.
In[46]:=
Click for copyable input
Out[46]=
定义一个函数来计算并且显示 , , ... 等等的前几个特征值，以及范数界限.
In[46]:=
Click for copyable input
In[48]:=
Click for copyable input
Out[47]=
范数界限在这个例子中非常明显.

Korteweg-deVries

Korteweg-deVries 偏微分方程是在浅水表面的波浪的数学模型：
我们考虑如下边界条件：
并且在区间 上求解.
采用直线法的离散化用来形成由 192 个常微分方程组成的方程组.

步长

选择要求解的问题.
In[49]:=
Click for copyable input
在 LSODA 中使用的向后差分公式方法在求解这一问题时遇到了困难.
In[50]:=
Click for copyable input
Out[50]=
一个图形显示了步长快速减小.
In[51]:=
Click for copyable input
Out[51]=
相反， 立即切换到使用线性隐式欧拉方法，它需要很少的积分步骤.
In[52]:=
Click for copyable input
Out[52]=
In[53]:=
Click for copyable input
Out[53]=
一旦在积分开始的时候选择了刚性方法，插值法从来不会切换回非刚性求解器.
因此，这是一个非刚性检测的最差情况例子的形式.
尽管如此，使用子空间迭代的成本仅是总积分时间的百分之几.
在关闭到非刚性方法的转换时，计算所采用的时间.
In[54]:=
Click for copyable input
Out[54]=

雅可比实例

采用前面定义的监控函数，收集前几个雅可比矩阵的数据.
In[55]:=
Click for copyable input
初始雅可比 的图形显示.
In[57]:=
Click for copyable input
Out[57]=
计算并且显示 , , ... 的前几个特征值，以及范数界限.
In[63]:=
Click for copyable input
Out[57]=
范数界限略被高估，但更重要的是，他们没有提供实部和虚部的相对大小的显示：

选项总结

StiffnessTest

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

方法选项 的选项.

NonstiffTest

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

方法选项 的选项.

Ask a question about this page  |  Suggest an improvement  |  Leave a message for the team
格式:   HTML  |  CDF