微分代数方程

微分代数方程（differential-algebraic equation，DAE）是同时包含微分关系和代数约束的方程组。它不仅描述状态怎样变化，也规定哪些状态根本不允许出现。电路中的电流守恒、刚性连杆的长度限制，都可能使若干变量不能自由指定。

状态演化与即时约束

考虑无量纲教学方程组 $x^{'} = - x + z, 0 = x + z - 1.$ 第二式要求任意时刻都满足 $x + z = 1$ ，所以 $z$ 不是可以独立指定初值的第二个自由状态。若给定 $x (0) = 0$ ，就必须有 $z (0) = 1$ ；若同时指定 $z (0) = 0$ ，则初值已经违反约束，无连续解能在零时刻满足全部方程。

消去代数变量 $z = 1 - x$ ，得到 $x^{'} = 1 - 2 x$ 。利用初值求解并恢复 $z$ ： $x (t) = \frac{1 - e^{- 2 t}}{2}, z (t) = \frac{1 + e^{- 2 t}}{2} .$ 两者逐渐靠近 $1 / 2$ ，而和始终为一。把这两个式子代回，既要检查微分式，也要检查代数式；只验证一个分量的导数还不够。

解在直线x加z等于一上从零一移动到二分之一二分之一 — 代数约束把可行状态限制在直线 x+z=1 上。实线段标出该初值的未来轨迹，箭头朝向平衡点；标出的不相容初值位于直线外。

隐式形式与可解性

一般隐式形式写作 $F (t, y, \dot{y}) = 0, y \in ℝ^{n} .$ 若 $F$ 连续可微，且关于导数的雅可比矩阵 $F_{\dot{y}}$ 可逆，隐函数定理可局部把它化成普通的显式微分方程 $\dot{y} = f (t, y)$ 。真正带代数约束的情形，通常存在不通过这一步解决的奇异结构；仅凭写成隐式等号，还不能断言问题困难。

常用的半显式形式是 $x^{'} = f (t, x, z), 0 = g (t, x, z) .$ 其中 $x \in ℝ^{m}$ 是微分变量， $z \in ℝ^{k}$ 是代数变量。若 $f, g$ 足够光滑，并且方阵 $g_{z}$ 在解附近可逆，隐函数定理给出 $z = ϕ (t, x)$ ，代回即可得到局部常微分方程。此时称为典型的指数一半显式系统。

对约束求一次时间导数，可以看出代数变量的导数被怎样确定： $0 = g_{t} + g_{x} x^{'} + g_{z} z^{'}, z^{'} = - g_{z}^{- 1} (g_{t} + g_{x} f) .$ 这一步不是把原约束删掉；原约束仍负责选择正确的解曲面。SUNDIALS 的 IDA 文档说明了通用隐式残量形式与相容初值的计算。^[1]

相容初值与隐含约束

开篇系统还有导数关系 $x^{'} + z^{'} = 0$ 。在初始点 $(0, 1)$ ，原微分式给出 $x^{'} (0) = 1$ ，进而得到 $z^{'} (0) = - 1$ 。因此初始状态与初始导数也不能任意组合。对更多约束的系统，即使位置满足约束，速度也可能不相容。

以二维平面内长度为 $ℓ > 0$ 、质量为一的理想摆为例。位置为 $(q_{1}, q_{2})$ ，竖直向上取正，速度为 $(v_{1}, v_{2})$ ， $g > 0$ 为重力加速度。用乘子 $λ$ 表示沿杆方向的约束力系数： ${q_{1}}^{'} = v_{1}, {q_{2}}^{'} = v_{2}, {v_{1}}^{'} = - λ q_{1}, {v_{2}}^{'} = - g - λ q_{2}, q_{1}^{2} + q_{2}^{2} = ℓ^{2} .$ 位置约束求导，得到速度必须切于圆： $q_{1} v_{1} + q_{2} v_{2} = 0.$ 再求导并代入加速度，得到 $v_{1}^{2} + v_{2}^{2} - g q_{2} - λ ℓ^{2} = 0, λ = \frac{v_{1}^{2} + v_{2}^{2} - g q_{2}}{ℓ^{2}} .$ 位置和速度一致后，乘子的初值才由此确定。若要进一步得到全部未知量、包括 $λ$ 的显式导数，还需再微分一次约束关系；按通常微分指数约定，这个笛卡尔坐标形式是指数三系统。把“第二次微分可算出乘子”误说成指数二，会混淆算出 $λ$ 与算出 $λ^{'}$ 。微分指数的定义及摆例可参见 Atkinson、Han 与 Stewart 的教材。^[2]

摆的杆长约束使位置在圆上，速度箭头沿圆切线 — 圆周是位置约束，切线箭头是满足 q·v=0 的瞬时速度。法向约束力负责使加速度满足杆长要求；速度并不沿杆方向。

保持约束的离散计算

对开篇例子，在新时刻同时满足隐式欧拉微分式和代数式： $\frac{x_{n + 1} - x_{n}}{h} = - x_{n + 1} + z_{n + 1}, x_{n + 1} + z_{n + 1} = 1.$ 消去 $z_{n + 1}$ 后得到 $x_{n + 1} = \frac{x_{n} + h}{1 + 2 h}, z_{n + 1} = 1 - x_{n + 1} .$ 取 $h = 1 / 4$ 、 $x_{0} = 0$ ，第一步为 $(x_{1}, z_{1}) = (1 / 6, 5 / 6)$ ，第二步为 $(5 / 18, 13 / 18)$ 。每步代数约束都精确成立，但这不等于时间解精确；在 $t = 1 / 2$ ，真实 $x = (1 - e^{- 1}) / 2 \approx 0.316060279$ ，数值值 $5 / 18 \approx 0.277777778$ 仍有时间离散误差。

一个变式是把约束改为 $x + z = 2$ ，保持 $x_{0} = 0$ 。此时消元为 $x^{'} = 2 - 2 x$ ，精确解 $x = 1 - e^{- 2 t}$ ，递推变为 $x_{n + 1} = (x_{n} + 2 h) / (1 + 2 h)$ 。同样的约束结构可以复用，右端常数和初始代数变量却必须一同改变。

一般 DAE 求解器会在每步求解非线性残量，使用雅可比矩阵，并控制误差；这里的显式消元只适用于这个小例子。在高指数问题中，直接微分约束后丢弃原约束会产生漂移：连续恒等式只规定约束导数为零，数值舍入和离散误差却可能让约束常数偏离原值。降指数、约束投影或使用合适的结构保持方法，需要与具体系统匹配。

参考资料

↑ SUNDIALS Development Team. IDA: Mathematical Considerations，§5.2.1 Initial Condition、§5.2.2 IVP solution。
↑ Kendall Atkinson, Weimin Han, David Stewart. Numerical Solution of Ordinary Differential Equations，第10章，第163—165页的摆和微分指数，§10.1（第165—168页）的相容初值与约束漂移。

[ida-1] SUNDIALS Development Team. IDA: Mathematical Considerations，§5.2.1 Initial Condition、§5.2.2 IVP solution。

[iowa-2] Kendall Atkinson, Weimin Han, David Stewart. Numerical Solution of Ordinary Differential Equations，第10章，第163—165页的摆和微分指数，§10.1（第165—168页）的相容初值与约束漂移。

[1]

[2]