一根刚杆既能运动,又必须保持长度。描述它的方程因此不仅规定速度怎样变化,还规定当前状态能出现在哪里。给每个变量随意填一个初值,可能根本没有轨道能从那里出发。
形式陈述
方程、解和初始数据
令 I ⊂ R 为开时间区间,x ( t ) ∈ R d 。隐式微分系统写成
(1) F ( t , x ( t ) , x ′ ( t ) ) = 0 , F : Ω ⊂ R × R d × R d → R d . 本页把它作为微分代数方程 的总框架,重点讨论关于导数变量的Jacobian 理路 Jacobian 矩阵 Jacobian matrix 多元映射各偏导数组成并表示其导数的矩阵。 F v 可能奇异的情形。若它在某个满足方程的三元组附近可逆,隐函数定理可将式(1)局部改写为显式ODE 理路 常微分方程 Ordinary differential equation · ODE 未知函数及其单一自变量导数组成的方程。 x ′ = a ( t , x ) 。
经典解 是逐点满足式(1)的 C 1 曲线。初始状态 x 0 称为一致 ,若存在包含 t 0 的时间邻域及经典解,满足 x ( t 0 ) = x 0 。若同时提交导数 v 0 ,还要求该解满足 x ′ ( t 0 ) = v 0 。
代入一次得到 F ( t 0 , x 0 , v 0 ) = 0 只是必要条件。约束在随后每个时刻仍要成立,对它求导可能产生额外的初始条件。相容性是一条轨道能否存在的要求,不只是某次残差计算为零。
可逆代数块给出的局部定理
一个重要形式是
(2) x ′ = f ( t , x , z ) , 0 = g ( t , x , z ) , x ∈ R n , z ∈ R m , n , m ≥ 1. 其中 x 称为微分变量,z 称为代数变量;这描述所给方程中的角色,不是变量永久不变的身份。
设 f , g 在 ( t 0 , x 0 , z 0 ) 的开邻域为 C 1 ,并满足
可 逆 (3) g ( t 0 , x 0 , z 0 ) = 0 , g z ( t 0 , x 0 , z 0 ) 可逆 . 则通过 ( x 0 , z 0 ) 存在唯一局部经典解,并且初始导数唯一确定为
(4) x 0 ′ = f ( t 0 , x 0 , z 0 ) , z 0 ′ = − g z − 1 ( g t + g x f ) ( t 0 , x 0 , z 0 ) . 唯一性限于所选点附近的同一解支。定理没有保证轨道永久存在,也没有允许它穿过 g z 退化的位置。
直觉
先解约束,再解运动
对 g ( t , x , z ) = 0 使用隐函数定理 理路 隐函数定理 Implicit function theorem 当相关偏导块可逆时,方程组局部可把部分变量表示为其余变量的函数。 ,把 ( t , x ) 一并当作参数,得到唯一 C 1 函数
z = ζ ( t , x ) , g ( t , x , ζ ( t , x ) ) = 0. 代回第一式,只需解
(5) x ′ = f ( t , x , ζ ( t , x ) ) . 右端为 C 1 ,因此对状态局部Lipschitz,Picard–Lindelöf定理 理路 Picard–Lindelöf 存在唯一性定理 Picard-Lindelof theorem · Cauchy-Lipschitz theorem 连续且对状态变量局部一致 Lipschitz 的向量场给出常微分方程初值问题的唯一局部解。 给出唯一局部 x 。再令 z = ζ ( t , x ) ,便恢复完整解。反过来,附近的任一DAE解都必须位于这张约束图上,因此也满足同一ODE,唯一性随之成立。
沿解对约束求导得到
0 = g t + g x x ′ + g z z ′ . 用 x ′ = f 解出 z ′ ,就是式(4)。即使原方程没有写 z ′ ,约束仍规定了它怎样随时间变化。
为什么不能只保留求导后的方程
如果只解
(6) x ′ = f ( t , x , z ) , z ′ = − g z − 1 ( g t + g x f ) , 链式法则仅保证 d g ( t , x , z ) / d t = 0 。这意味着约束残差保持初值,不意味着它等于零。初始残差为 c ≠ 0 时,式(6)会沿着错误的水平集 g = c 前进。
一致初始化把轨道放到正确的约束面;积分过程中重新检查原约束,则防止数值误差悄悄改变这个水平集。这是两项不同的工作。
例子与边界
一个非线性约束和两个代数分支
考虑
(7) x ′ = − x + z , z 2 − x − t = 0 , ( x ( 0 ) , z ( 0 ) ) = ( 1 , 1 ) . 这里 g z = 2 z = 2 ,所以附近选择正分支 z = x + t 。初始数据满足约束,且
x ′ ( 0 ) = 0 , z ′ ( 0 ) = 1 2 . 若把同一个 x ( 0 ) = 1 配成 z ( 0 ) = − 1 ,也得到一致状态,但属于另一局部分支,此时 x ′ ( 0 ) = − 2 、z ′ ( 0 ) = 1 / 2 。可逆偏导保证的是每个给定状态附近唯一,并非整个代数方程只有一个根。
即时残差为零,仍可能没有经典解
取
(8) x ′ = z , x = t 2 . 在 t 0 = 0 提交 x 0 = 0 , z 0 = 1 , x 0 ′ = 1 ,两条即时方程都成立。但对 x = t 2 求导要求 z = 2 t ,所以任何经典解都必须有 z 0 = 0 。提交状态不一致。
这里约束对代数变量的偏导 g z = 0 ,式(3)不适用。进一步求导会产生哪些条件,由微分指标与隐藏约束 理路 微分指标与隐藏约束 Differentiation index of a DAE · DAE differentiation index · 隐藏约束 · 微分指标 以导数阵列提取全部状态导数所需的最少微分次数定义指标,逐层推导指数1、2及机械指数3的隐藏条件,并展示约束降阶后的漂移和数据微分放大。 具体分析。
奇异导数块本身不是困难程度的结论
方程 ( x ′ − 1 ) 3 = 0 的导数Jacobian在解上为零,却与简单ODE x ′ = 1 完全等价。这是退化的方程表示,不是实际增加了三层约束。
另一个边界是 z 2 = x 在 ( 0 , 0 ) 处:可逆代数块条件失败,附近有两支;不能把只在正分支内证明的结论跨过该点使用。一般奇异系统还可能无解、欠定或发生秩变化,式(1)本身并不附带存在唯一性保证。
推论与应用
后向Euler怎样同时满足两组方程
沿用隐式时间步 理路 隐式 ODE 方法与非线性步求解 Implicit ODE method · Implicit time stepping · Implicit time integration 把隐式时间步分解为离散方程、非线性求解与线性求解,区分时间稳定性、迭代收敛、阻尼及误差预算。 ,从一致状态 ( x n , z n ) 出发,在 t + = t n + h 求
(9) R 1 = X − x n − h f ( t + , X , Z ) = 0 , R 2 = g ( t + , X , Z ) = 0. 不必先显式求出 ζ ,也不能仅更新 X 而把旧 z n 当成新约束值。Newton矩阵为
( I − h f x − h f z g x g z ) . 若 g z 可逆,用Schur消元 理路 Schur 补与静态凝聚 Schur complement · Static condensation 把内部变量的精确响应折入界面矩阵与右端,推导 Schur 补、恢复公式及最小能量性质,并逐项凝聚五节点链。 消去代数修正,微分块恰为
(10) S = I − h ( f x − f z g z − 1 g x ) . 括号正是约化右端 f ( t , x , ζ ( t , x ) ) 对 x 的导数。因此足够小的步长下,邻近解支存在唯一;也可直接在 h = 0 对式(9)用隐函数定理,此时块三角Jacobian的对角为 I , g z 。
非线性迭代并非因此从任意初猜都收敛。实现仍须保留邻近分支、检查 g z 条件性,并分别记录两块残差。稠密求解可分解 m 阶代数块及 n 阶Schur块;成本还含多右端求解和矩阵乘法,直接组装的一轮为 O ( ( n + m ) 3 ) 量级,函数及导数评值另计。
把一个离散步完整复算
对式(7)取 h = 1 / 4 ,式(9)成为
5 X = 4 + Z , Z 2 = X + 1 4 . 消元得到 20 Z 2 − 4 Z − 21 = 0 。与 Z = 1 连续相接的根为
(11) Z = 1 + 106 10 , X = 41 + 106 50 . 另一根落在负分支,不是本初值的小步延续。将式(11)代回两组残差,它们都精确为零;这认证了后向Euler步,不代表已经得到原DAE在 t = 1 / 4 的精确轨道值。
大型模型还可能给出奇异质量矩阵 E x ′ = A x + f 。此时正则常系数DAE 理路 正则常系数微分代数系统 Regular linear DAE · Constant-coefficient differential-algebraic system · Weierstrass decomposition for DAEs 从可逆移位与核像分解证明正则矩阵束的微分—幂零分块,给出全部经典解和一致初值,并把外力导数损失与隐式步的误差放大写成有限证书。 会区分能够自由设定的微分初值和被外力及其导数强制确定的代数分量;单凭某一步矩阵可逆,仍不足以判断输入初值一致。
参考资料