引言

一般来说，常微分方程组（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 符号式求解.

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

利用选项 "MassMatrix" 生成形如 的系统：

利用选项 "Residual" 生成形如 的系统：

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

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

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

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

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

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

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

    .

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

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

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

微分代数方程的解法

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

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

• IDA -（SUNDIALS 程序包中的隐式微分代数求解器），基于根据反微分公式（BDF）.

• StateSpace - 隐式求解 的导数 ，以使用底层的常微分方程求解器.

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

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

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

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

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

具有不变量的常微分方程

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

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

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

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

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

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

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

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

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

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

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

    .

利用选项 "MassMatrix" 求解形如 的系统：

微分代数方程的指标约简

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

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

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

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

        钟摆力的示意图.

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

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

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

指标约简方法.

下列选项可用于所有指标约简方法.

指标约简选项.

指标约简通过对微分代数方程组的方程进行微分完成. 假定在指标约简的过程中对方程 eqn 进行微分，使得 deqn=eqn 包括在系统中 eqn 的位置. 一旦微分完成，微分后的方程将组成一个完全确定的系统，并能够自行求解. 然而，在通常情况下非常重要的是将原始方程作为约束纳入系统中. 这是由 "ConstraintMethod" 指标约简选项控制的.

约束方法.

比较两个解对约束 的符合程度：

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

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

指标约简算法

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

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

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

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

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

结构式关联

获取结构式关联.

获得方程组 、 对于变量 和 的结构式关联矩阵：

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

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

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

获得方程组 、 对于变量 和 的结构式矩阵：

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

获得方程组 、 对于变量 和 的结构式矩阵，其中对第二个方程进行了微分：

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

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

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

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

这是一个全矩阵，因为 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 算法的不同及其优越性，考虑线性系统

    .

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

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

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

    .

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

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

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

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

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

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

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

约束方法

虚拟导数

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

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

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

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

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

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

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

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

以 为变量的关联矩阵是

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

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

试着对系统进行积分，将手动引入的虚拟导入替换 和 . (使用虚拟导数是默认设置，因此不必给出选项设置.)

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

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

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

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

投影

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

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

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

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

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

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

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

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

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

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

    .

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

    .

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

    .

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

    .

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

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

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

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

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

    ,

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

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

    .

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

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

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

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

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

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

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

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

配点法

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

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

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

    ,

    .

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

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

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

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

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

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

"Collocation" 方法选项.

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

配点方向

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

"Forward"	配点在范围 内进行
"Backward"	配点在范围 内进行
"Centered"	配点在范围 内进行

"CollocationDirection" 子选项.

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

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

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

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

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

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

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

配置点

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

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

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

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

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

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

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

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

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

配置范围

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

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

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

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

"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 算法最多需要处理一个大小为 的矩阵，而不是配点法的大小为 的矩阵. 在大多数情况下，该算法只需要处理远小于 的矩阵. 这使得该方法在处理非常大的系统时，相当有效. 本节将提供底层算法的简要说明.

已知一个形如 的指标约简系统，可以首先将系统关于起始猜测 和 线性化为

.

在雅可比矩阵 上执行 QR 分解，使得 ，其中 是一个正交矩阵， 是上三角矩阵， 是一个置换矩阵. 对于微分代数方程，矩阵 将不会有满行秩. 根据矩阵 的秩，将线性化系统写作

.

对变换后的系统分两步求解：

.

一旦找到分量，将使用置换矩阵，把变换后的变量转换回原来的形式，并且不断重复该过程，直到收敛.

QR 方法可以从方法 "DAEInitialization"，利用Method->{"DAEInitialization"->{"QR",qropts}} 调用，其中 qropts 是 QR 方法的选项.

QR 初始化的选项.

考虑单摆的例子.

由于这是一个高指标的微分代数方程，必须对系统进行指标约简. 一旦系统指标被约简为1/0，必须对所有的变量及其导数进行初始化.

在前面的例子中，给出的初始条件给会导致一组一致的初始值. 在初始条件不一致的情况下，该算法试图找到一致的初始值.

通过指定选项 "DefaultStartingValue"，迭代所用的起始猜测被改变了. 对于不一致的初始条件，这将导致不同的一致集合.

BLT 方法

BLT 方法基于的思想是，将指标约简系统变换到一个分块下三角形式，然后迭代式求解系统的子集. 子集自身使用高斯-牛顿方法求解. 这有时是有利的，因为生成指标约简系统的指标约简算法往往已经有了一个分块下三角形式，而对于较小的子系统，相应的矩阵计算也较小，通过处理这种子系统，减少了计算的复杂性. 该算法还可以检查不一致的初始条件，并试图对其纠正.

BLT 方法可以利用 Method->{"DAEInitialization"->{"BLT",bltopts}}，从方法 "DAEInitialization" 调用，其中 bltopts 是 BLT 方法的选项.

BLT 方法的选项.

要了解 BLT 方法的优点，重要的是要认识到，指标约简算法导致一个扩展的指标为0或1的方程组. 考虑下述化学系统作为例子，该系统描述了一种放热反应.

原来的微分代数方程由四个方程组成，而指标约简后的系统包括八个方程. 指标约简后的方程如下.

注意，扩展后的系统现在包含虚拟导数. 对扩展的系统进一步分析，可以看到系统的结构.

这表明，扩展的系统实际上可以按顺序来求解. 因此，可以按顺序求解单个方程，而不是求解由九个方程组成的方程组. 从初始化的角度来看，设置 t->0 将允许变量的初始化.

这是利用 BLT 方法进行初始化的主要原理. BLT 方法处理具有下三角形式的大型系统时，尤其有用，并且计算效率高.

对于初始条件不一致的情况，BLT 方法可以适当修改初始条件. 但是，需要注意的是，由 BLT 方法和 QR 方法获得的结果，可能会有所不同. 考虑单摆具有不一致的初始条件的例子.

不一致的初始条件

微分代数方程系统组成及其求解所面临的主要挑战之一是，起始的初始条件必须是一致的. 一开始，这看起来似乎很简单，但对于高指标的微分代数方程系统，很容易遇到隐藏的条件/方程，初始条件必须同时使其成立. 如果 NDSolve 遇到它认为不一致的的初始条件，会试图纠正这些问题.

为了说明隐藏约束和不一致初始条件这一点，考虑摆系统(2).

首先考虑第一种情况，即系统初始条件违反代数约束时. 这可以通过故意设置初始条件为 ，因此它总是违反约束 .

如果将该系统提供给 NDSolve，将产生 NDSolve::ivres 的信息，表明初始条件被认为是不一致的. 然后，它试图找到一个一致的初始条件，使用该不一致的值作为其迭代的起始猜测值.

有时，可能并不是很清楚指定条件一致与否. 有可能初始条件看似满足约束了，但可能无法满足隐藏约束.

这些初始条件满足方程 ，但不满足隐式约束 . 这种约束并不明显，因为它是隐藏的，并且只能通过指标约简发现. 配点算法不执行指标约简，并作用于原方程。即时是在指标约简系统的情况下，配点算法也能检测到初始条件的不一致，并尝试纠正它.

具有事件的微分代数方程重新初始化

考虑一个由两个方程组成的系统，该系统含有一个离散变量，其中状态在 处发生变化.

在前面的示例中， 和 的初始条件已知. 初始化代码自动计算出一致的初始条件 . 在时刻 ，事件被触发，离散变量应该从1变化到2. 在时刻 ，需要进行变量的重新初始化. 变量 被称作修正变量（modified variable）. 重新初始化过程采取以下策略：

1.  保持所有的状态变量、离散变量和修正变量固定. 如果状态变量具有固定的导数（或指定要更改），那么该变量可以被改变.

2.  如果策略1重新初始化失败，则保持状态变量和修正变量固定，并允许其他变量及其导数改变.

3.  如果策略2失败，则修正变量和微分代数方程中的具有导数的状态变量保持不变. 所有其他变量及其导数可以被改变.

4.  如果策略3失败，那么只有修正变量保持不变，允许所有的状态变量及其导数改变.

重新初始化之前，、 和 的值如下.

获取事件发生之前，在 的值：

基于上述策略概述，在时刻 ，使用策略1和2，它要求、 和 不改变. 由于 ，这显然违反代数方程 . 因此，策略1和2失败.

策略3要求具有导数的变量不变. 因此 保持不变（同时允许 ) 和 改变. 这改变了 的值，如下图所示（蓝色线条）.

获取事件发生之后，在 的重新初始化值：

WhenEvent 允许对事件行为进行同时和顺序计算，包括状态变化.

在每个状态变化后，即使对于相同的事件，也会使用前策略对微分代数方程重新初始化，所以如果可能的话，只给出一条状态变化所需的规则更为有效. 最终状态可能有所不同，因为有多个状态变化在同一事件时间完成，而每一个后续计算和设置将使用新的值，对于微分代数方程，这些值将取决于重新初始化. 考虑前面的问题，其中附加了其他事件.

得到的解与前面的不同. 如要了解原因，必须按照顺序设置. 首先将 设为2，从而在计算其他设置之前，对微分代数方程重新初始化. 重新初始化将 的值变成了2，原因是该策略这样能够得到与余项一致的解. 执行第二个行为，将 设置为 2（ 的当前值），但这与余项不一致，因此，由于 是一个代数变量，它的值被设置为4 (). 最后，设置 为 ，将 重置为2，为了避免修改新设置的代数变量， 的值被重置为1.

查看变化顺序的一种方法是，在两个设置之间添加 Print 命令.

因此，不同的设置顺序将导致不同的解.

对于单一规则，只存在一种初始化. 因此，如果设置了所有变量，所赋的值应该是一致的，否则重新初始化策略将失败，或者不能给出期望的值.

这里，, 和 是固定的. 这是一个非常严格的约束，所有的策略都失败了.

将一个变量设置为它的当前值，在一个单一规则内有效，强制它保持不变.

在 时，将 变为2，强制 保持不变：

由于 保持不变， 被强制改变了.

微分代数方程范例

圆环问题

这是一个取自 Campbell 和 Moore [CM95]的检验问题. 微分代数方程描述了三维粒子的运动. 这是一个指标为3的微分代数方程，具有三个位置变量、三个速度变量和一个约束. 解的形式是四维的，确切的解由下式给出：

    .

从油箱中排出压力

这是石化工程中动态处理的一个简化模型. 关于该模型的详细信息，请参阅 Hairer et. al. [HLR89]. 这是一个指标为2的系统.

刚性曲柄滑块

这个例子是一个简单的约束多体系统. 该系统由一个质量为 m1、铰接在固定表面上的曲柄、一个质量为 m3 的滑块和一个质量为 m2 的连接杆组成. 滑块被限制为仅在水平面内移动. 由于各个组件质量的存在，在系统构建过程中必须考虑惯性矩. 该问题的示意图如下. 该系统可使用两个角度完全定义.

曲柄滑块机制的状态完全可以由两个角度 和滑块与原点的距离 完全描述.

通过求解力并应用牛顿定律 和 ，推导运动方程. 曲柄和连杆具有质量，因此也有惯性.

定义系统的物理参数. 这里， 和 是惯性矩.

弹簧阻尼器附加在链接上的曲柄滑块

这是一个简单曲柄滑块机构的延伸. 在曲柄和连杆之间安装了一个弹簧和阻尼器系统. 此弹簧约束着曲柄的运动范围，而阻尼器则缓慢地使系统达到稳定状态. 系统示意图如下.

曲柄滑块机制的状态可以完全由两个角度 和滑块距原点的距离 描述.

通过求解力并应用牛顿定律 和 ，推导运动方程. 曲柄和连杆具有质量，因此也有惯性.

双连杆平面机械臂

双连杆平面机械臂是一个具有约束的无摩擦机械系统. 该系统由两个链接组成，这两个链接通过一个关节彼此相连，使得第一个连杆的一端保持固定，而第二个连杆的端部的轨迹是指定的. 该系统的示意图如下所示.

该系统可以完全由链接之间的角度定义. 关于该系统的更多信息，请参阅 Ascher et. al. [ACPR94].

定义机械臂端部必须遵循的轨迹. 这对应于定义 的位置. 运动为

通过求解力并应用牛顿定律推导运动方程. 关于角度的微分方程也可推导得到. 因此，笛卡尔坐标 和 也以角度的形式表示.

n 摆问题

构造一个 摆系统的函数：

晶体管放大电路

阿克苏诺贝尔化学反应

下面的系统描述了一种化学过程，此过程中有两个化学品 FLB 和 ZHU 在 不断添加到系统的情况下进行混合. 化学名称都是虚构的. 系统的详细信息以及与它相关联的动力学信息可以从[N93]得到. 反应方程式如下. 反应速率是指定的，它与每一个反应相关联. 代数约束来自最后的反应，这需要 FLB、ZHU 及其混合物保持恒定的关系.

汽车轴问题

汽车轴模型是一个简单的多体系统，其目的是模拟在颠簸的道路上行驶的汽车轴线的行为. 系统的示意图如下所示. 该模型是由两个轮子组成：左轮被假定为在平坦的路面上移动，而右轮则移到了不平坦的路面. 崎岖不平的表面用正弦曲线建模. 通过两个无质量的弹簧将左右车轮的运动传输到汽车底盘，如下图所示. 汽车底盘用一个质量为 M 的长条模拟. 关于该模型的更多信息，请参阅 Mazzia et. al. [MM08]. 这是一个指标为3的系统.

左轮胎总是在平坦的表面上，而右轮胎则周期性地经过一个凸块.

车轮之间的距离必须始终保持固定，并且汽车的底盘必须始终为一个固定的长度.

一维组合式椭圆抛物线偏微分方程

下面是一个合成示例，求解以非线性方式耦合在一起的抛物线和椭圆偏微分方程. 偏微分方程的空间离散化导致一个含有显著数量个代数变量的系统. 方程的设置使得解析解存在于系统. 以下给出耦合系统的详细信息.

多物种的食物网模型（二维组合式椭圆抛物线偏微分方程）

这个例子是一个多物种的食物网模型 [B86]. 与典型的捕食-食饵问题不同，这个问题有一个空间分量. 捕食-食饵关系在二维空间域上通过一系列扩散参数定义. 以下给出一般的耦合偏微分方程：

    .

在这个系统中，前 个物种为食饵，而接下来的 个物种为捕食者. 这里考虑 的具体情形. 这表示有一个诱饵和一个捕食者. 空间域用20×20的网格进行离散化. 这意味着最终的微分代数方程组有2×400个变量. 最初指定捕食者的分布为均匀分布. 这导致一个不一致的初始条件，这可以由内置算法纠正.

在 和 处对解进行可视化：

安德鲁斯的榨汁机理

安德鲁斯的榨汁机理是一个多体系统，包含7个不同的链接. 各个链接的尺寸及其规格如下. 该系统可使用七个角度完全定义. 角度在下面的系统中标记为 q₁–q₇. 此外，有六个与链接尺寸和连接相关联的代数约束. 这导致了一个维度为14、指标为3的系统. 该系统的详细信息可以参阅 Hairer 和 Wanner [HW96].

与各个链接相关联的维度如下图所示. 七个链接被标记为 K₁–K₇.

用于此模型的参数记号及其说明可参阅 Hairer 和 Wanner [HW96].

系统的运动可以用七个角度来定义. 这里，角度定义为 . 运动方程可以利用哈密顿原理推导，得到一个形如 的方程组，其中 是一个质量矩阵， 是力的矢量， 是约束的矢量，λ 是拉格朗日乘子.

定义与系统相关联的质量矩阵 ：

定义与系统相关联的各种力 ：

定义与各种链接相关联的约束矢量 ：

建立方程 和具有特定参数值的约束方程：