微分代数方程的数值解
引言
一般来说，一个常微分方程组（ODEs）可以表示为正规形式，
应变量 
的导数使用自变量 
和应变量 
显式表示. 只要函数 
具有足够的连续性，我们就总是可以找到一个初值问题的解，其中对自变量的一个特定值，给出应变量的值.
对于微分代数方程，导数一般没有显式表示. 事实上，一些应变量的导数通常没有在这些方程中出现. 微分代数方程组的一般形式为
其中关于 
的雅可比 
可能是奇异的.
微分代数方程组可以通过对自变量 
求导来转化为常微分方程组. 一个微分代数方程的指数（index）实际是需要对微分代数方程组求导以获得常微分方程组的求导次数. 尽管求导是可能的，它一般不作为一个计算技术，因为原有的微分代数方程的属性常常在求导后的方程的数值模拟中丢失.
因此，微分代数方程的数值方法被设计来对微分代数方程组的一般形式使用. NDSolve 中的方法一般设计来求解指数-1 微分代数方程组，但是可能对于高指数问题也适用.
本教程将用大量的例子来展示求解微分代数方程和常微分方程的一些不同之处.
微分代数方程的初始条件的指定与常微分方程的大不相同. 对于常微分方程，正如我们已经讨论过的，一组初始条件唯一地确定了一个解. 对于微分代数方程，情况远非这么简单；找到满足方程的初始条件甚至都可能很困难. 为了更好地理解这个问题，考虑下面的例子 [AP98].
以下是具有三个方程的微分代数方程组，但只具有一个微分项.
初始条件明显不是自由选取的的；第二个方程要求 
或者是 0 或者是 1.
以下从
的导数的一个特定的初始条件开始求解微分代数方程组.
| Out[12]= |  |
为了得到这个解，NDSolve 首先使用 Solve 和与 FindRoot 很相似的程序的组合来搜索满足这些方程的初始条件. 一旦找到相符的初始条件，微分代数方程就可以使用 IDA 方法求解.
| Out[13]= |  |
以下显示解的图形. 解
被解
遮蔽了，它具有相同的常量值1.
| Out[15]= |  |
不过，对于所有满足方程初始条件，解也有可能不存在.
以下尝试从具有导数0的稳定状态开始，使用
寻找一个解.
| Out[16]= |  |
| Out[17]= |  |
如果查看把 
设置为 1 时的方程，用户就可以看到为什么不可能在 
后继续进行求解.
把
代入方程.
| Out[18]= |  |
中间的方程实际上消失了. 如果用户在 
时对最后一个方程求导，用户获得条件 
，但这时是第一个方程与初始条件中 
的值不一致.
结果，
的唯一解是 
，由该解可知，系统具有指数 2.
该问题的另一组解是当 
时. 用户可以通过把这指定为一个初始条件来求解.
以下使用
寻找一个解. 我们也必须对
指定一个值，因为它是一个微分变量.
| Out[19]= |  |
| Out[21]= |  |
一般说来，用户必须对微分变量指定初始条件，因为通常这里有一个含参数的通解（parametrized general solution）. 对于满足 
的该问题，通解是 
，所以有必要给出 
来确定解.
NDSolve 不总是可以找到与方程一致的初始条件，因为有时候这是一个困难的问题. "通常，应用中求解一个微分代数方程组最难的部分是确定一组相容初始条件来开始计算". [BCP89]
| Out[22]= |  |
如果 NDSolve 未能找到相容的初始条件，用户可以使用具有一个良好的初值的FindRoot 或者一些其它的程序来获得相容的初始条件. 如果用户已知接近一个良好的起开始猜想值的值，则 NDSolve 使用这些值开始搜索，这可能会有所帮助. 用户可以指定应变量的值以及它们的导数.
使用微分代数方程组的指数-1 系统，常常可以先求导并且使用一个常微分方程求解器来获得解.
以下是罗伯逊化学动力学问题. 由于大的和小的变化率常数，这个问题相当刚性.
以下通过对平衡方程求导，把罗伯逊动力学问题作为常微分方程求解.
| Out[26]= |  |
和 
的主要变化具有两个完全不同的时间度量， 决定了该问题的刚性.
以下显示解
和
.
| Out[27]= |  |
| Out[33]= |  |
一个给定变量的解看上去将很接近，但是比较化学平衡约束表明了它们之间的不同.
以下是平衡方程的常微分方程解和微分代数方程解的误差图形，分别用黑色和蓝色显示. 这里使用一个 log-log 标度，因为
以及误差的变化都很大.
| Out[37]= |  |
这种情况下，在预期的容差以外，两个解都很好的满足平衡方程. 请注意，虽然平衡方程的微分代数方程解的误差在某些点比较大，从长期来看，一旦通过了快速变化范围，微分代数方程解就会回转，从而更好地满足约束.
用户可能想要求解形式如下的一些微分代数方程
即微分方程的解需要满足一个特定的约束. NDSolve 不能直接处理这样的微分代数方程，因为指数太高，而 NDSolve 期望方程数目与应变量数目相同. 不过，NDSolve 具有一个通常可以求解该问题的 Projection 方法.
这样的一个约束系统的一个很简单的例子是模拟单摆运动的非线性振子.
以下定义方程、不变量约束和起始条件来模拟单摆的运动.
注意，微分方程实际是不变量的导数，所以求解该方程的一个方法是使用不变量.
以下使用不变量方程求解一个单摆的运动.
选项告诉
NDSolve 不要符号地求解
的二次方程，而是把该系统作为一个微分代数方程求解.
| Out[58]= |  |
不过，这可能不是我们所期待的解：当从 
开始时，不变量方程具有 
常数的解. 事实上，从这点开始并没有唯一的解. 这是因为如果你去求解 
，函数将不满足唯一性所要求的连续性.
以下只使用微分方程求解单摆的运动方程. 这里使用了方法
，因为它也可以是投影(projection)方法的子方法.
| Out[59]= |  |
| Out[60]= |  |
| Out[62]= |  |
不变量的误差不是很大，但是它确实显示了一个稳定一致的漂移. 最后，它可以大到影响解的保真度.
以下求解单摆的运动时，把每个步骤中的运动限制在不变量上.
| Out[63]= |  |
以下显示
NDSolve 使用投影法，在所采用的时间步的末尾不变量的图形.
| Out[65]= |  |
