微分代数方程的数值解—Wolfram 语言参考资料

link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

引言

一般来说，常微分方程组（ODE）可以表示为正规形式：

因变量的导数以独立瞬态变量和因变量的形式显式表示. 只要函数具有足够的连续性，总能找到初值问题的唯一解，其中对自变量的一个特定值，给出因变量的值.

对于微分代数方程，导数一般没有显式表示. 事实上，一些因变量的导数通常不出现在这些方程中. 例如，下列方程

不显式包含的任何导数. 这种变量通常被称为代数变量.

微分代数方程组的一般形式为

求解微分代数方程组往往需要许多步. 下面的流程图说明了与 NDSolve 中求解微分代数方程相关的一般步骤.

在 NDSolve 中求解微分代数方程组所涉及的步骤流程图.

一般来说，微分代数方程组可以通过对自变量求导，转换为常微分方程组求解. 微分代数方程的指标（ index ）是需要对微分代数方程组求导以获得常微分方程组的求导次数.

内置于 NDSolve 中的微分代数方程求解方法对指标为1的方程组有效，对于指标较高的方程组则需要约简指标以得到解. 可以指示 NDSolve 如何约简指标. 如果发现方程组的指标高于1， NDSolve 将生成一则信息及求解微分代数方程所需的一组步骤.

求解高指标微分代数方程的第一步是约简方程组的指标. 通过求导约简指标的这个过程叫做指标约简（index reduction），可以在 NDSolve 中用符号式过程完成.

指标约简的过程产生一个新的等价方程组. 在其中一个方法中，这个新的方程组是通过将微分变量的位置用新的（虚拟）变量替代而重新组成的. 这产生一个能够被唯一求解的扩展系统. 另一种重新建立方程组的方法涉及到对原始方程组求导次，直到所有变量的微分得到处理. 前面的个方程组被视为不变量.

为了求解新的指标约简方程组，必须找到一组一致的初始条件. 规范形式的常微分方程组总能通过给出初始时刻的值来初始化. 然而，对于微分代数方程，得到满足残余方程的初始条件可能会比较困难；这相当于求解非线性代数方程组，其中只有一些变量可以独立地指定. 这将在后面的章节中详细讨论. 此外，初始化必须是一致的. 这意味着残余方程的导数也需要成立. 通常，方程组的指标越高，初始化越困难. 在这种情况下， NDSolve 一般不了解的组分如何作用，因此不能进行自动的指标约简. 但是， NDSolve 提供有一定数目的不同方法，可以通过选项访问，进行指标的约简，并得到微分代数方程的一致性初始化.

本文档其余部分的结构方式如下：首先介绍与微分代数方程指标相关的一些术语. 随后的一节处理指标降一次的微分代数方程的求解. 然后介绍更高阶指数的微分代数方程的指数约简方法. 随后讲解微分代数方程的一致性初始化. 最后一节将在前面各节给出的范例的基础上，给出更多的范例.

这里载入范例将会用到的程序包：

微分代数方程的指标

对于一般的微分代数方程，系统的指标指的是以和的形式求解，系统中部分或全部方程所需的最小求导次数. 注意文献中有多种关于指标的不同概念；本文将使用的是上述定义.

考虑下列微分代数方程

为了得到关于的微分方程，必须对第一个方程进行一次求导. 这将得到新的系统

求导后的系统现在含有必须考虑的 . 要得到关于的方程，第二个方程必须求导三次. 这将得到

由于部分方程必须求导三次才能得到关于的导数公式，该微分代数方程的指标为3. 最终得到的求导系统被称为指标为0的系统. 通过进行适当的代换，可以查看通过求导得到的常微分方程，在本例情况下，从第一个中减去第二个方程. 将得到的方程可以写为

同样地，指标为1的系统可以通过对第二个方程两次微分，得到以下系统

这是典型的，将系统的指标降为1是非常常见的做法，因为底层积分例程可以有效地处理它们.

微分代数方程的指标可能并不总是很明显. 考虑这样一个系统，根据初始条件，该系统可能会出现指标为1或2的行为.

有三个方程组成的微分代数方程组，但只有一个微分项：

对于这个微分代数方程的初始条件不是自由的；第二个方程要求为0或1.

如果，则其他两个方程构成一个指标为1的系统，因为对第三个方程微分得到的是两个常微分方程的系统

另一方面，如果，则其他两个方程构成一个指标为2的系统. 对第三个方程微分得到，代入第一个方程得到，这需要被微分以得到常微分方程

由于且，初始条件中也有没有额外的自由度.

求解该微分代数方程组，首先指定

：

所需的初始条件数和微分代数方程的指标均会对所得到的实际解产生影响.

微分方程的处理

为了求解微分方程， NDSolve 将用户指定的系统转换成的三种形式之一. 这一步是很重要的，因为系统的构造方式不同，所选择的微分方法也会不同.

通过使用选项 Method -> { "EquationSimplification"-> simplification } ，可以指定方程的表示形式. 下面的化简方法，可以由 "EquationSimplification" 选项指定.

Automatic 尝试以

的形式符号式求解系统. 如果不能求解或求解时间太长，将使用 "MassMatrix" 或 "Residual" 方法化简. 这是默认设置. "Solve" 尽可能以

的形式符号式求解系统. "MassMatrix" 尽可能将系统化简为

的形式. "Residual" 从方程左端减去方程右端，形成残差函数

方程化简方法.

考虑系统

在默认情况下，将尝试以自变量和因变量的形式显式求解微分 . 出于提高效率的目的， NDSolve 将首先利用 LinearSolve 试着求解显式系统. 但如果失败的话，系统将利用 Solve 符号式求解.

NDSolve 会自动尝试将系统显式表示：

当符号解能够得到时，使用 "Solve" 方法与默认方法的效果是相同的.

利用选项 "MassMatrix" 生成形如

的系统：

利用选项 "Residual" 生成形如

的系统：

当系统为残差形式时，导数由一组新的变量表示，这些产生的变量符号是唯一的. 这些工作变量和指定变量之间的对应可以利用 NDSolveStateData 的属性 "Variables" 和 "WorkingVariables" 得到.

显示工作变量和指定变量之间的对应关系：

产生显式常微分方程组的过程，有时会由于系统的符号处理而变得很费时。由于这个原因，获得方程有一个1秒的时间约束. 如果超出这个时间， NDSolve 会生成一条信息，并尝试将系统作为微分代数方程求解.

定义非线性方程组的变量和方程：

利用默认的方程化简策略求解：

下列方法选项可用来指定化简方法，以便更好控制 "EquationSimplification" 的行为.

选项名称

默认值

"TimeConstraint"1 允许 Solve 显式求导的最长时间（秒） "SimplifySystem" False 是否化简，通过为尽可能多的因变量获取解析解

正如这条信息所示，由于 Solve 超出了默认的1秒时间约束， NDSolve 未能对导数显式求解. NDSolve 获得显式表达式所需的时间可以利用选项 "TimeConstraint" 控制.

利用子选项 "TimeConstraint" 控制 NDSolve 获得显式表达式所需的时间：

在上述范例中，由于第一个方程的平方项，我们得到了四个可能的解. 当作为微分代数方程求解时，数值解将根据其初始化方式选择其中的一个解的分支.

为了演示方法选项 "SimplifySystem" 的用法，请尝试下面的例子：

上面方程组的分析解可以通过把代入第一个方程，求得的解，然后代入第二个方程组，给出的解. 方程组的最后解为

在没有首先执行方程组的指标化简以及形成新的等价常微分方程组前， NDSolve 不能求解该方程组. 这个代价是昂贵的. 但是当开启子选项 "SimplifySystem" ， NDSolve 可以检测到可以求得的解析表达式/解的变量，并执行重复替代回原方程组. 这将导致原始方程组被化简或在某些情况下（如上所示）返回解析解.

用 "SimplifySystem"-> True 求解上面的微分代数方程：

如果对于某些变量，不存在解析解， NDSolve 则会返回插值函数来作为解.

求解一个方程组，其中，只有对某些因变量可以求得解析解：

开启子选项后，常量参数也可以直接代回方程组.

求解包含常量参数的常微分方程组：

微分代数方程的解法

在 NDSolve 中内置有多种求解微分代数方程的方法.

两种方法可以处理指标为1的微分代数方程的一般残差形式：

IDA -（ SUNDIALS 程序包中的隐式微分代数求解器），基于根据反微分公式（BDF）.
StateSpace - 隐式求解的导数，以使用底层的常微分方程求解器.

如果系统的指标高于1，上述两种求解器通常都会失败. 如果精确求解这类系统，需要进行指标降低. 注意到指标为1的系统可以降低为常微分方程，但使用上述求解器往往更为有效.

NDSolve 也有一些求解器，可以求解能够转换为特殊形式的微分代数方程.

MassMatrix - 求解形如的微分代数方程

Projection - 求解具有不变量、形如的常微分方程

下面的小节将介绍使用这些方法求解微分代数方程中的一些方面.

具有不变量的常微分方程

经常会遇到下列形式的微分代数方程

其中是一个与微分方程一致的不变量.

当微分代数方程组通过指标降低转换为常微分方程的形式时，被微分以得到常微分方程形式的方程与常微分方程是一致的. 并且对解进行积分时，确保这些方程完全成立通常非常重要.

这类系统具有的方程数比未知量多，因而是超定的，除非这些约束与常微分方程一致. NDSolve 不能直接处理这里微分代数方程，因为它需要系统的方程数与因变量个数相同. 为了求解这类系统，内置于 NDSolve 的 "Projection" 方法对不变量进行处理，在每一个时间步后将计算得到的解投影到 . 这确保了代数方程随着解的演化而始终成立.

一个关于这种约束系统的例子是模拟钟摆运动的非线性振荡器.

这里定义方程、不变量约束以及一个模拟钟摆运动的初始条件：

需要注意的是的微分方程实际上是不变量的导数，所以解方程的一个方式是使用这个不变量.

如返回的信息所示， NDSolve 不直接求解该系统：

使用 "Projection" 时间积分法得到了在各时间步均满足不变量的解：

绘制各时间步的不变量；它被满足的精度接近于机器精度：

求解具有质量矩阵的方程组

当导数表现为线性形式时，微分系统可以写成以下形式：

其中是一个常被称作质量矩阵的矩阵. 有时，也被称为阻尼矩阵. 如果矩阵是非奇异的，则可以被转换，从而使系统作为常微分方程求解. 然而，如果奇异矩阵存在，系统可以作为微分代数方程求解. 在任何情形下，利用这种特殊形式是有利的.

表示圆柱体上质点运动的系统

可以用下式以质量矩阵表示

定义系统变量、方程和初始条件：

利用选项 "MassMatrix" 求解形如

的系统：

微分代数方程的指标约简

NDSolve 中微分代数方程的求解器能够完全自动地处理指标为 1 的微分代数方程组. 对于指标更高的系统，则有必要进行指标降低以得到解. 指标降低可以通过给以 NDSolve 适当的选项来进行.

NDSolve 使用符号方法进行指标降低. 这意味着如果你的微分代数方程组的形式为，其中 NDSolve 不能看到的分量如何交互和计算符号导数，指标约简不能完成，并且因为这个原因，在默认情况下 NDSolve 不能进行指标约简.

以描述钟摆运动的经典微分代数方程为例：

其中，在点处有一个质点，被长度为的绳子约束. 是拉格朗日乘子，实际上就是绳子的拉力. 为简单起见，在指标约简的介绍中，取 . 下图显示的单摆系统的原理图.

钟摆力的示意图.

如果系统的指标在方程中不明显，我们需要在 NDSolve 中尝试一下，如果该求解器不能求解，在许多情况下，它能够产生一条消息，指示在指定的初始条件下表现为什么指标.

为约束单摆系统定义方程以及一组完整的初始值：

试着求解约束单摆系统：

这则信息表示指标为 3，意味着指标约简对于求解系统是必要的. 注意这里使用了完全一致的初始化，以避免可能产生的初始化问题. 关于求解步骤方面的内容在 “ 微分代数方程的一致初始化 ” 一节中有详细介绍.

选项设置 Method -> { "IndexReduction"-> { irmeth , iropts } } 用于指定 NDSolve 使用具有一般指标约简选项 iropts 的指标约简方法 irmeth .

None 不保留原始方程作为约束，仅求解微分方程 "DummyDerivatives" 使用 Mattson 和 Soderlind 方法 [ MS93 ] 将某些导数用代数变量（虚拟导数）替换 "Projection" 将指标约简为0以得到常微分方程，并使用 "Projection" 时间积分法，其中对原始方程微分作为不变量

约束方法.

使用指标约简并求解微分方程：

比较两个解对约束

的符合程度：

如果解只通过被微分方程得到，则摆的长度偏离1. 将原始方程作为约束纳入解是指标约简的一个重要方面.

下面两节将介绍指标约简算法，然后介绍可用于 NDSolve 的约束方法.

指标约简算法

NDSolve 具有两种进行指标约简的算法，Pantelides 算法和结构性矩阵算法. 这两种方法都是有方程的符号形式确定要微分的方程，然后使用符号微分得到微分后的系统.

两种方法都利用了方程中结构式关联的概念；特别地，如果变量显式出现在方程中，则称在中关联. 结构式关联矩阵是一个这样的矩阵：如果在中关联，其元素为；其他情况下 . 如果结构式关联矩阵可以重新排序，使得沿对角线的元素为1，则称微分代数方程组具有结构式指标1. 这种重新排序的存在意味着方程和变量是匹配的，即每个变量都存在一个方程.

Pantelides 算法适用于基于关联的图，对于极其庞大的系统这可以非常高效.

注意结构式指标可能会小于实际的系统指标. 在实数系统中，结构式存在的项可能会沿一些解消失，或者雅可比矩阵是奇异的. 结构矩阵方法试图考虑到可能的雅可比矩阵的奇异性和指数减少一些系统会做的更好的工作，但在更大的计算费用的成本. 结构式矩阵方法尽量考虑雅可比矩阵的奇异性，对于某些系统的指标约简效果较好，但需要较大的计算成本.

下一节将介绍获得结构式关联矩阵的命令，这对于接下来的一节所介绍的 Pantelides 和结构式矩阵算法是非常有用的.

结构式关联

NDSolve`StructuralIncidenceArray [ eqns , vars ] 返回 SparseArray ，对于方程 eqns 中的变量 vars 具有结构式关联矩阵的模式

获得方程组

的结构式关联矩阵：

注意 SparseArray 仅含有使用较少内存的非零模式，表明它代表关联模式. 它可以利用 ArrayRules 被转换为由1组成的 SparseArray .

定义一个命令，将模式数值转换为1和0：

系统具有结构式指标1，由于不需要重新排序以避免零在对角线上.

如果相对于导变量的结构式关联经处理可以避免零在对角线上，则系统的指标为0.

获得方程组

、

对于变量

和

的结构式矩阵：

没有重新排序，以避免对角线上的零. 但是，如果对第二个方程进行了微分，则是可能的.

获得方程组

、

对于变量

和

的结构式矩阵，其中对第二个方程进行了微分：

微分后的系统具有结构式指标0，证明原系统的指标为1.

考虑约束摆( 2) 作为另一个例子.

获取约束摆系统相对于所出现的最高阶导数的结构式关联矩阵：

矩阵不能重新排序以避免零在对角线上，因此指标大于1.

注意函数 NDSolve`StructuralIncidenceArray 只是对表达式的 FullForm 进行字面上的匹配，以检查关联. 它不对方程或变量的构成是否正确无误进行检查.

这里检验符号 x 和 y 在每个方程中的字面形式：

这是一个全矩阵，因为 x ' [ t ] 的 FullForm 是含有 x 的 Derivative [ 1 ] [ x ] [ t ] .

Pantelides 方法

Pantelides [ P88 ]提出的方法是一个图理论方法，最初提出的目的是为了得到微分代数方程一致初始化. 它适用于因变量和方程的二分图，当算法可以找到的图的遍历时，系统的指标即将到了最高为1. 在这个意义上，遍历实际上指的是对变量和方程排序，使图的关联矩阵的对角线上没有零.

当二分图没有完整的遍历时，该算法通过对方程进行微分而有效地扩充路径，在此过程中引入新的变量（先前变量的导数）. 除非原来的系统在结构上是奇异的，该算法将终止于遍历. 一般情况下，该算法仅为了得到遍历而进行必要的微分，但因为它是一种贪婪算法，它的数量并不总是最低的.

详细的算法描述不在本教程的范围之内. 内置于 Wolfram 语言相当密切地服从[ P88 ]介绍的算法. 它采用 Graph 有效地实现图计算，并在需要微分时用 D 对方程进行微分.

当系统遍历存在时，则有可能得到变量和方程的排序，使得关联矩阵为分块下三角（BLT）形式. BLT 形式用于设置虚拟导数法，以保持约束，也用于 "StateSpace" 时间积分法中. 在 "微分代数方程的状态空间方法" 中对 BLT 排序进行了说明.

在约束摆系统( 2) 中开始指标约简，首先考虑变量 . 一般要从出现在方程的最高阶导数开始. 在此提醒您，摆方程是

相对于最高阶导数的图关联矩阵为：

关联矩阵的构建通过检查变量在方程中的存在实现. 如果变量存在于方程中，则权重为1；否则为0. 因此，对第一个方程中的变量进行检查，得到（正如在关联矩阵的第一行所看到的）.

由第一个关联矩阵可以发现遍历不存在. 对最后一个方程微分两次

现在，关联矩阵为

这可以重新排序，通过交换第一个和最后一个方程，使得对角线上没有零，从而得到关联矩阵

所得的的系统的指标为1，因为仍作为代数变量出现，并且在微分后的系统中不含有导数. 对系统再次进行微分，将为提供一个微分系统. 因此，原系统的指标确实为3.

作为第一种方法，经过两次微分的方程可以直接求解.

求解具有两次微分方程的系统：

绘制应为不变的量，基于最初指定的方程：

原摆长的约束条件没有满足，事实上，它伸长了. 没有满足约束条件的原因是新的系统没有等价地表示原来的系统. 这可以通过新系统不含的原约束，而是它的微分后的方程这一事实看出来.

为了用新系统模拟原来的动态系统，额外的方程必须得到满足. 然而，包括额外的方程会导致方程数多于未知数个数. 为了解决这个问题，可使用虚拟导数的方法，以使得指标约简的系统变得平衡.

Pantelides 算法是有效的，因为它使用的图形算法即使对于非常大的问题也具有良好控制的复杂度. 然而，由于它是基于纯粹的关联，可能出现导致雅可比矩阵奇异的问题. 结构式矩阵方法可能解决这种问题，但对于大型系统可能会具有较大的计算复杂度.

结构式矩阵算法是 Pantelides 算法的替代方法. 结构式矩阵法遵循从 Unger [UKM95] 和 Chowdhary 等 [CKL04] 人的工作. Pantelides 算法等基于图的算法依靠遍历来约简指标. 然而，这种算法不考虑具体系统中的变量有时可能消去的情况. 这导致算法低估了系统的指标. 如果微分代数方程没有正确地约简到指标0或指标1系统，那么数值积分可能会失败或产生不正确的结果.

作为 Pantelides 方法没有考虑变量消去的例子，考虑下面这个微分代数方程

要进行指标约简，作为第一步，对第二个方程求导一次，得到：

Pantelides 算法在第一次求导后停止，原因是它得到了变量的导数和遍历. 由于所有导数都是由第一次求导得到，Pantelides 算法估计系统指标为0. 指标为0的系统等价于常微分方程组. 求导后系统的导数的雅可比矩阵为

显然是奇异的. 通过简单替换，有可能看到求导后的系统等价于

这意味着为了得到常微分方程组，必须对第二个方程进行第二次求导. 这就得到了最终的系统

上述系统对于导数具有非奇异的雅可比矩阵. 现在上述系统具有正确的指标0，由此可以得到一个常微分方程组. 由于需要两次求导，系统的实际指标为2（而不是1）.

与 Pantelides 算法不同，结构式矩阵方法考虑与线性项关联的所有消去. 该方法的步骤利用作为下述一阶系统的单摆范例 ( 2) 进行说明：

第一步构建两个矩阵和 . 两个矩阵分别是导数与和变量相关的关联矩阵. 对于上述系统，给出矩阵如下：

上述矩阵含有一些元素. 这里用表示占位符，表示该方程（行）的变量在方程中以非线性形式出现. 注意到第二个方程，项和作为乘积出现. 因此，第二行的第一列和最后一列用标记.

一旦矩阵构建好了，下一个目标是将微分代数方程组转换为常微分方程组. 为此，矩阵必须是满秩的. 为实现这个目的，首先确定需要求导的方程. 在这种情况下，最后一个方程需要求导. 在求导后，结构式矩阵为

由于求导，矩的最后一列改变了. 为更清楚起见，考虑 ( 3) 中最后一个方程. 对方程求导给出 . 注意到项在方程中出现，并且是以非线性形式. 因此，与导数相关的关联行是，与变量相关的关联行是 .

下一步涉及到对矩阵以高斯消去的形式进行矩阵分解. 的最后一行中的元素，可以利用行1和3消去. 这种消去也影响矩阵的行. 将的第一行乘以，并减去最后一行，得到

将第三行乘以 *，并从最后一行减去，得到

矩阵的秩仍然不足，因此必须重复上述步骤. 第二次迭代后得到系统

第三次迭代后得到

矩阵现在是满秩矩阵. 因此迭代停止. 由于需要对进行三次迭代才得到满秩矩阵，该系统的指标为3. 值得注意的是，与 Pantelides 方法不同，结构式矩阵算法作用于关联矩阵.

下一个例子将展示结构式矩阵方法与 Pantelides 算法的不同及其优越性，考虑线性系统

设置方程和初始条件：

试着利用 Pantelides 算法对系统求解：

失败的原因是存在一个奇异的雅可比矩阵. 当微分代数方程的指标大于1时，求解器发出信息 NDSolve ::nderr 或 NDSolve ::ndcf 是非常常见的.

利用结构式矩阵方法进行指标约简，从而求解系统：

精确的解是及 . 比较精确结果和数值结果，可以发现计算得到的解是相当准确的.

接下来，考虑对系统指标过低或高估计的情况. 考虑这个线性方程组：

设置的方程、初始条件和变量。

所有变量的精确解是 . 这是一个指标为4的系统. Pantelides 算法所估计的指标为2.

试着利用 Pantelides 算法进行指标约简，以求解系统：

该求解器遇到了收敛失败，原因是雅可比矩阵是奇异的.

利用结构式矩阵算法进行指标约简，以求解系统：

NDSolve 可执行某些检验，确定要用的指标约简方法. 对于这种情况，使用选项 Automatic ， NDSolve 选取的是结构式矩阵算法.

设置 "IndexReduction"-> Automatic ，以对系统求解：

如果系统被强制视为指标0，则可以执行积分. 然而，在指标估计不正确的情况下，结果可能会明显地偏离，因为解是通过一种完全不同的方法计算的.

利用 Pantelides 算法将系统指标当作零进行求解：

由于该问题的精确解是已知的，希望看一下数值解和精确解之间的差异.

由于对系统的指标估计过低，没有进行必要的求导，从而导致对解的不准确逼近. 相反地，结构式矩阵方法能够正确确定系统的指标.

作出由 "StructuralMatrix" 得到的系统的解的图形：

利用结构式矩阵方法，将系统指标作为0进行求解：

结构式矩阵方法确实能更好地处理这种系统. 但是需要注意的是，这种考虑消去的做法的严谨性是要在其他方面付出代价的. 结构式矩阵法依赖生成、维护和操作矩阵，与 Pantelides 算法相比，这需要占用显著的计算开销. 因此，这种算法仅适用于中小规模的问题.

约束方法

虚拟导数

虚拟导数法的目的师引入代数变量，从而使得原方程组和为了指标约简而求导后的方程所组合的系统不会超定.

再次考虑单摆系统( 2) .

这里利用将虚拟导数内置于 NDSolve 的方法进行指标约简，从而对系统求解：

注意在这个例子中包括了 "ConstraintMethod"->"DummyDerivatives" 是为了强调，而不是必要的，因为尝试用虚拟变量是默认设置.

根据原来指定的方程绘制应为不变的量：

现在约束条件被很好地满足了.

Mattson 和 S ö derlind [ MS93 ] 的虚拟导数法的主要思想是通过引入表示导数的新变量，使得再次引入约束方程而不会得到超定系统成为可能.

回想一下系统( 2) 的指标约简方程是

其中有两个约束来自第三个方程和它的第一个导数，

令和分别为替换和的代数变量. 由于包括了约束，这些变量的整个系统有下列各式给出：

以为变量的关联矩阵是

这个矩阵可以重新排序，使得对角线上的元素非零.

对关联矩阵重新排序：

一般而言，该算法可以应用于任何指标约简的微分代数方程；然而，在对系统进行积分时，求解器可能会遇到奇异性的情况.

试着对系统进行积分，将手动引入的虚拟导入替换

和

. (使用虚拟导数是默认设置，因此不必给出选项设置.)

上述系统不能完成积分，因为系统在时变成了奇异的.

这里将它所能计算的解进行绘图：

从图形可以看出， NDSolve 在变为 0 时被卡住了. 当为 0 时，多个项消去了，关联矩阵成为奇异的，也不能再找到非零的对角线，因此虚拟导数系统实际上在沿着解通过时指标 .

考虑另一种虚拟导数替换方法，即用和替换和 . 遗憾的是，具有这种替换的系统在处是奇异的.

补救的方法是使用动态选择，当远离0时，用和替换和，当变得过小时，用和替换和，动态切换至系统. 在必要时， NDSolve 自动使用具有状态规则的 WhenEvent 来处理具有虚拟导数的动态状态选择. 动态状态选择的精确细节超出了本教程的范畴. 详尽的解释可以在 Mattson 和 Soderlind [ MS93 ] 中找到.

投影

利用虚拟变量的替代方法是将指标约简为 0，得到一个常微分方程组，并使用投影来确保原约束被满足.

对于单摆系统 ( 2)，系统的指标为 3. 这意味着为了获取关于所有变量的常微分方程，必须对第一个和第二个方程求导一次，而必须对第三个方程求导三次. 这将得到

能够求解得到导数、和 .（ NDSolve 将自动将系统约简至第一阶.）

为了正确表示系统的动态，必须考虑第一个求导后的方程. 这些方程是

上述方程被视为必须在每个时间步成立的不变方程. 由于常微分方程组是通过对这些方程求导得到的，从而确保约束是与常微分方程一致的，因此利用 "Projection" 时间积分法，将约束视作不变量，能够对系统求解.

利用 "Projection" 方法求解常微分方程：

这可以直接在 NDSolve 中利用内置的指标约简方法完成.

这里利用内置于 NDSolve 的投影进行指标约简，从而求解系统：

通过这种方法，原约束通常能够在不超过 NDSolve 的局部限度内被满足，这由 PrecisionGoal 和 AccuracyGoal 选项指定.

系统具有的额外不变量也应该被满足. 其中一个是能量守恒定律，可表示为 . 由于投影方法只在时间步长的端部进行投影，这是做比较最合适的地方.

这里得到在时间步的解的能量，并对虚拟变量的能量和投影解进行绘图：

偏微分代数方程（PDAE）的指标约简：

考虑下面的偏微分方程组：

在上面的方程组中，第二个 PDE 不包含任何变量的时间导数分量. 因此， NDSolve 不能使用行方法求解该方程组，因为由此导出的质量矩阵会是奇异矩阵并且方程组具有指标1. 为了减少方程组的指标，上面方程组要表示为矩阵向量形式：

表示离散化空间导数得到的矩阵，表示在离散空间区间，离散化变量得到的向量. 指标约简在新的方程组中使用上述的指标约简方法进行. 由此导出的指标约简方程组为：

现在，空间导数引回方程组：

导出的方程组可以用 NDSolve 中传统的方法求解.

定义方程组：

对 PDAE 进行指标约简求解方程组：

注意：对于空间离散化使用了200点，因为基于常数初始化条件的默认空间网格不足以处理随时间发展的解的变化性.

绘制解：

值得注意的是，指标约简是在时间导数下进行的，因此无需或添加新的边界条件. 然而，如果 PDAE 是高指标的，那么就需要额外的初始化条件.

微分代数方程的一致初始化

微分代数方程组(DAEs)的最一般表示形式为

其中可包括微分方程和代数约束. 为了得到的解，需要一组和的一致初始条件，以开始积分. 一致性的一个必要条件是初始条件满足然而，这个条件自身可能并不充分，原因是对原方程求导所生成的新的方程也需要满足初始条件. 因此这个任务就是试着找到一组满足所有一致性条件的初始条件. 寻找一致初始条件可以说是求解微分代数方程的最困难的部分之一.

考虑下面的线性微分代数方程问题：

微分代数方程含有的导数. 这将表明，为了唯一求解这个问题，必须指定的初始条件. 然而，您不能给变量指定任意的初始条件，因为存在一组唯一的一致初始条件. 对该代数方程求导一次，得到代入第二个方程，得到的解. 对求导，并代入第一个方程，得到的解. 这表明微分代数方程组的解是固定的，因此初始条件也是固定的，并且的仅有的一组一致初始条件是重要的是要注意，关于初始条件的约束不是很明显，一般使得找到一致初始条件这个问题相当具有挑战性.

在某些情况下，您可能只是想要得到微分代数方程问题的一致初始化，而不是进一步计算解. 这可以用 NDSolve 通过指定时间积分的端点与初始条件所在时间是相同的实现. 当您只需要检查开始的条件，而无需执行完整的积分时，这样的步骤一般是有帮助的.

这里计算微分代数方程组的初始条件：

NDSolve 提供了微分代数方程初始化的多种方法. 方法 m 可以用 Method -> { "DAEInitialization"-> m } 指定.

Automatic 自动确定要使用的方法 "Collocation" 使用配点法（Collocation）；该算法可用于高指标系统的初始化 "QR" 对于指标为1的系统，使用基于 QR 分解的算法 "BLT" 对于指标为1的系统，使用分块下三角排序（BLT）法

微分代数方程初始化的方法.

配点法专门用来处理微分代数方程组，将其作为残值黑盒，并试图在初始点附近的一小段区间上满足残值. 此功能使该方法可用于处理高指标系统的初始化. 然而，由于该算法不分析微分代数方程组，该系统的结构不能被利用以提高计算效率. 这使得该算法比其他方法慢.

QR 和 BLT 方法专门用于求解指标为1和指标为0的系统. QR 法依靠分解导数的雅可比矩阵生成两个解耦的子系统，从而使得变量及其导数的初始条件可以迭代计算. 这使得该方法非常有效和强大.

BLT 法依赖于检查理微分代数方程组的结构，把原来的系统分割成许多较小的子系统，从而有效地求解每一个子系. 这种方法要求处理小得多的雅可比矩阵（与整个系统相比），从而对于大型系统的计算是有效的.

以下各节详细介绍 NDSolve 如何利用不同的算法，得到微分代数方程的一致初始条件.

配点法

基本的配点法将因变量以基函数的形式展开为

为得到和，试着强制条件在时间为的个配置点成立. 为了实现这个目的，隐式微分代数方程被线性化为

其中和表示分别相对于和的雅可比矩阵：

将线性化应用于个并置的时间点，得到系数的线性方程组，这可以用牛顿法迭代求解.

其中通过直线搜索法得到. 如果有个系数和个并置点，则您将处理的是一个 ( 5 ) 中的 × 的线性方程组.

对于没有进行指标约简的高指标系统， NDSolve 将自动使用配点算法来初始化系统. 考虑单摆系统( 2)

为了开始对系统积分，必须知道的初始条件. 对于这种系统，将表示一些物理性能，可以指定为合理值. 然而，的值可能不是那么容易找到. NDSolve 可以在这里使用，以确定适当的值.

获取摆的微分代数方程的一组一致的初始条件：

请注意，只给出了部分初始条件，该算法正确地计算所有所需的初始条件.

配点法的默认设置已被选择，处理各种不同系统的初始化，但是，在某些情况下，更改设置是必要的，以找到一个一致的初始化. 下列选项可用于调整算法.

配点法选项名称

默认值

"CollocationDirection" Automatic 选择配置点所要放置的方向；可能的选项包括 "Forward" 、 "Backward" 、 "Centered" 、 "Automatic" "CollocationPoints" Automatic 所用的配置点数；使用子选项 "ExtraCollocationPoints" ，将添加更多的网格点，而同时保持固定的逼近阶 "CollocationRange" Automatic 配置点分布的范围 "DefaultStartingValue" Automatic 用于未指定的组件中的非线性迭代的初始值；选项可以是单元素值 ig / { ig } 或数对 { ig 1 , ig 2 } 的列表，其中 ig 1 , ig 2 分别是变量及其导数的开始值 "MaxIterations"100 要执行的最大迭代次数

"Collocation" 方法选项.

接下来的章节将更详细地讨论这些选项，并显示在什么样的情况下，它们可能会帮助找到一致的初始条件.

配点方向

为了进行初始化，选择一组配置点. 通过选项 "Centered" 、 "Forward" 或 "Backward" ，点可以分别放置在感兴趣的点的两侧、前面或后面. 下表描述的是选择这些选项之一的效果.

"Centered" 配点在范围

内进行

"CollocationDirection" 子选项.

上面的指的是计算初始条件的时间. 是偏移值. 缺省情况下， NDSolve 自动选择配点方向.

配点方向在处理不连续性时是非常重要的.

在这里，在时刻，不连续.

NDSolve 自动设置 "CollocationDirection" 的选项.

应用 "Centered" 于 "CollocationDirection" ，将导致解发散.

然而，利用 "Forward" 或 "Backward" 方向，可以避免不连续的导数，从而得到正确的收敛.

在大多数情况下， NDSolve 自动确定应该采用的方向. 但是，如果已经知道在一个特定的方向上有可能存在不连续，则可以通过此设置避免. 反过来，这有助于减少 NDSolve 的计算时间.

配置点

选项 "CollocationPoints" 指的是配置点的个数，正如( 5 )中所示，用于级数逼近( 4) . ( 4) 中级数的阶数根据配置点来确定. 级数的阶数为 . 配置点在本质上是分布在配点范围上的切比雪夫点.

为了说明 "CollocationPoints" 的效果，考虑下面这个指标为3的例子.

先前的微分代数方程的解可以解析地找到. 精确一致的初始条件如下所示.

"CollocationPoints" 方法找到正确初始化的默认设置：

手动设定配置点是可能的. 设置 "CollocationPoints"->2 等效于( 4) 中的设置 . 这意味着该解使用常数和线性基函数逼近得到.

使用常数和线性基函数进行初始化：

尽管迭代收敛，结果也不是正确的一致初始条件. 因为它们充分满足指定的方程，因此不会发出任何消息.

检验指定方程中的初始条件：

还记得对于高指数方程，一致的初始条件应该也适用于求导后的系统.

增加配置点数，以允许更高阶的基函数，从而使算法收敛至期望的结果.

试着求解一系列的配置点的初始化问题：

将结果作为配置点的函数进行绘图，指示解如何收敛. Red 、 Green 和 Blue 分别对应变量 x 1 、 x 2 、 x 3 ：

由于逼近的阶数由配置点确定，( 4) 所得到的线性化的系统将导致一个方阵，并且由此得到的线性方程组可以有效地使用 LinearSolve 求解. 然而，产生的线性系统很可能是病态或奇异的. 这样的系统是敏感的，因此，可导致不稳定的迭代. 要解决这个问题，往往需要使用 LeastSquares 求解超定方程组.

使用子选项 "ExtraCollocationPoints" ，级数逼近的顺序保持与 ( 4) 中的相等，但 ( 5 ) 中的值修正为，其中是额外的配置点. 因此，这就产生了 ( 5 ) 中的超定方程组，可使用最小二乘法求解.

初始化一个超定系统：

配置范围

对于指标非常高而没有经过指标约简的系统，可能需要使用高得多的配置阶数，以及一个规定的范围，在该范围内进行配置. 缺省情况下，根据扩展所用的工作精度和基函数数目计算配置范围.

这里用一个指标为6的系统的作为范例，这个系统可以找到一个精确解.

定义指标为6的线性微分代数方程组的变量和方程：

在这种情况下，利用 DSolve 得到精确的解：

制作一个初始条件的表格：

初始化的默认选项对某些项给出不正确的结果：

通常情况下，使用一个高于系统指标的阶数进行逼近，通常是必要的. 由于这是一个指标为6的系统，必须依赖于 ( 4) 中的高阶逼近. 正如上一小节中提到的那样，这可以利用选项 "CollocationPoints" 做到.

利用15阶逼近的足够的点进行初始化：

解找到了，但与准确一致的初始条件相比有显著误差. 造成这种大误差的原因是，在这个例子中，默认的点之间距离太近，所以有明显的数值舍入误差.

将15个配置点遍布于长度为1的配置范围，进行初始化，从而避免数值的舍入误差：

"CollocationPoints" 和 "CollocationRange" 之间有密切的关系. 配置范围根据配置点的数目计算，从而避免过度的舍入误差，满足较高精度的计算. 除非明确指定， "CollocationRange" 计算为，其中是 WorkingPrecision 选项的设置，是 "CollocationPoints" 方法选项的设置.

指定初始条件

常微分方程和微分代数方程之间的主要区别之一是，为了对常微分方程组进行积分，必须提供所有变量的初始条件. 另一方面，对于微分代数方程，这些初始条件中的一些可能是固定的，并且必须直接从系统计算. 为了详细说明，考虑下面的常微分方程

前述系统的初始条件可以任意指定. 现在考虑相同系统的一个变形

该系统是一个微分代数方程，但是与常微分方程组具有相同的动态. 现在的主要区别是，初始条件和不是固定的，而是由代数方程确定.

前面的例子足够简单，可以确定哪些变量是固定的，也许需要手动计算初始条件. 然而，随着系统变得更加复杂和庞大，就必须依靠额外的工具.

在某些情况下，初始条件完全由方程组决定. 为了说明这种情况，考虑下面的线性微分代数方程组

虽然该系统包含两个导数，解和初始条件完全由系统决定. 不需要更多的信息.

初始化系统，不指定初始条件：

得到一致的初始条件：

值得注意的是，往往不可能知道哪些条件是固定的，哪些条件是自由的. 洞察正确初始条件或初始条件阶数的一种方法是，直接把没有初始条件的微分代数方程放在 NDSolve 中，处理 "DAEInitialization" 方法的选项.

考虑表示为一阶系统的单摆系统( 2) . 为了唯一地定义系统的状态，只有两个初始条件是必要的.

首先，考虑初始条件没有为微分代数方程给定的情况. 默认情况下，首先猜测所有没被指定的变量为1. 由于没有进行指标约简， NDSolve 自动选择配点算法.

可以观察到，该算法能够从无穷多组中找到一组可能的初始条件. 在这一点时，它当然是对起始猜测非常敏感的. 同样重要的是要注意，这不保证迭代的收敛性.

起始猜测可以使用配点法的选项 "DefaultStartingValue" 改变. 如果用 -1 作为起始猜测，则将得到一组不同的有效解.

正如前面提到的，两个初始条件足以固定系统的状态. 对于摆的情况，和可以固定系统.

计算摆系统的初始条件，其中摆的两个状态固定：

如果状态确实是固定的，则理论上讲，系统不应该对初始猜测敏感.

计算不同起始值的初始条件：

对于某些系统，初始条件的计算对起始猜测是敏感的. 为了说明这一点，考虑下面的系统.

对于这个问题 x ₂ [ t ]  0 ，通解是 { x ₁ [ t ]  x ₁ [ 0 ] + t ² / 2 , x ₂ [ t ]  0 , x ₃ [ t ]  t } ，因此有必要给出 x ₁ [ 0 ] 以确定解. 通过使用起始猜测值 1 ，强制算法改变指定的初始条件.

代码本质上切换到一个不同的求解途径.

出现这种情况的原因是，状态的默认起始猜测为 1 . 在牛顿迭代期间，由于第二个方程同样成立，系统变成奇异的，从而造成迭代发散. 为了避免完全失败，该算法被迫修改的指定条件.

为了避免这种行为，可以提供一个 “ 更好 ” 的猜测. 使用起始猜测 -1 ，迭代收敛到预期的一致初始条件.

QR 方法

QR 分解方法基于 Shampine 的工作 [S02] . 该算法专门用于求解指标为1或指标为0的系统. 该算法利用高效的 QR 分解技术对系统变量和导数进行解耦. 在这个过程中，每个迭代步骤的一个系统的解被细分成两个较小的方程组. 对于一个变量的微分代数方程组，QR 算法最多需要处理一个大小为的矩阵，而不是配点法的大小为的矩阵. 在大多数情况下，该算法只需要处理远小于的矩阵. 这使得该方法在处理非常大的系统时，相当有效. 本节将提供底层算法的简要说明.