Skip to content

模型Model

微分代数方程与一致初值

Differential-algebraic equation · DAE · 微分代数系统

将隐式微分关系、代数约束和一致初值放在同一模型中,以可逆代数偏导块证明局部可解性,并验证后向Euler的邻近分支和约束残差。

一根刚杆既能运动,又必须保持长度。描述它的方程因此不仅规定速度怎样变化,还规定当前状态能出现在哪里。给每个变量随意填一个初值,可能根本没有轨道能从那里出发。

形式陈述 ​

方程、解和初始数据 ​

令 I⊂R为开时间区间,x(t)∈Rd。隐式微分系统写成

(1)F(t,x(t),x′(t))=0,F:Ω⊂R×Rd×Rd→Rd.

本页把它作为微分代数方程的总框架,重点讨论关于导数变量的Jacobian Fv可能奇异的情形。若它在某个满足方程的三元组附近可逆,隐函数定理可将式(1)局部改写为显式ODE x′=a(t,x)。

经典解是逐点满足式(1)的 C1曲线。初始状态 x0称为一致,若存在包含 t0的时间邻域及经典解,满足 x(t0)=x0。若同时提交导数 v0,还要求该解满足 x′(t0)=v0。

代入一次得到 F(t0,x0,v0)=0只是必要条件。约束在随后每个时刻仍要成立,对它求导可能产生额外的初始条件。相容性是一条轨道能否存在的要求,不只是某次残差计算为零。

可逆代数块给出的局部定理 ​

一个重要形式是

(2)x′=f(t,x,z),0=g(t,x,z),x∈Rn, z∈Rm,n,m≥1.

其中 x称为微分变量,z称为代数变量;这描述所给方程中的角色,不是变量永久不变的身份。

设 f,g在 (t0,x0,z0)的开邻域为 C1,并满足

(3)g(t0,x0,z0)=0,gz(t0,x0,z0)可逆.

则通过 (x0,z0)存在唯一局部经典解,并且初始导数唯一确定为

(4)x0′=f(t0,x0,z0),z0′=−gz−1(gt+gxf)(t0,x0,z0).

唯一性限于所选点附近的同一解支。定理没有保证轨道永久存在,也没有允许它穿过 gz退化的位置。

直觉

先解约束,再解运动 ​

对 g(t,x,z)=0使用隐函数定理,把 (t,x)一并当作参数,得到唯一 C1函数

z=ζ(t,x),g(t,x,ζ(t,x))=0.

代回第一式,只需解

(5)x′=f(t,x,ζ(t,x)).

右端为 C1,因此对状态局部Lipschitz,Picard–Lindelöf定理给出唯一局部 x。再令 z=ζ(t,x),便恢复完整解。反过来,附近的任一DAE解都必须位于这张约束图上,因此也满足同一ODE,唯一性随之成立。

沿解对约束求导得到

0=gt+gxx′+gzz′.

用 x′=f解出 z′,就是式(4)。即使原方程没有写 z′,约束仍规定了它怎样随时间变化。

为什么不能只保留求导后的方程 ​

如果只解

(6)x′=f(t,x,z),z′=−gz−1(gt+gxf),

链式法则仅保证 dg(t,x,z)/dt=0。这意味着约束残差保持初值,不意味着它等于零。初始残差为 c≠0时,式(6)会沿着错误的水平集 g=c前进。

一致初始化把轨道放到正确的约束面;积分过程中重新检查原约束,则防止数值误差悄悄改变这个水平集。这是两项不同的工作。

例子与边界

一个非线性约束和两个代数分支 ​

考虑

(7)x′=−x+z,z2−x−t=0,(x(0),z(0))=(1,1).

这里 gz=2z=2,所以附近选择正分支 z=x+t。初始数据满足约束,且

x′(0)=0,z′(0)=12.

若把同一个 x(0)=1配成 z(0)=−1,也得到一致状态,但属于另一局部分支,此时 x′(0)=−2、z′(0)=1/2。可逆偏导保证的是每个给定状态附近唯一,并非整个代数方程只有一个根。

即时残差为零,仍可能没有经典解 ​

取

(8)x′=z,x=t2.

在 t0=0提交 x0=0,z0=1,x0′=1,两条即时方程都成立。但对 x=t2求导要求 z=2t,所以任何经典解都必须有 z0=0。提交状态不一致。

这里约束对代数变量的偏导 gz=0,式(3)不适用。进一步求导会产生哪些条件,由微分指标与隐藏约束具体分析。

奇异导数块本身不是困难程度的结论 ​

方程 (x′−1)3=0的导数Jacobian在解上为零,却与简单ODE x′=1完全等价。这是退化的方程表示,不是实际增加了三层约束。

另一个边界是 z2=x在 (0,0)处:可逆代数块条件失败,附近有两支;不能把只在正分支内证明的结论跨过该点使用。一般奇异系统还可能无解、欠定或发生秩变化,式(1)本身并不附带存在唯一性保证。

推论与应用

后向Euler怎样同时满足两组方程 ​

沿用隐式时间步,从一致状态 (xn,zn)出发,在 t+=tn+h求

(9)R1=X−xn−hf(t+,X,Z)=0,R2=g(t+,X,Z)=0.

不必先显式求出 ζ,也不能仅更新 X而把旧 zn当成新约束值。Newton矩阵为

(I−hfx−hfzgxgz).

若 gz可逆,用Schur消元消去代数修正,微分块恰为

(10)S=I−h(fx−fzgz−1gx).

括号正是约化右端 f(t,x,ζ(t,x))对 x的导数。因此足够小的步长下,邻近解支存在唯一;也可直接在 h=0对式(9)用隐函数定理,此时块三角Jacobian的对角为 I,gz。

非线性迭代并非因此从任意初猜都收敛。实现仍须保留邻近分支、检查 gz条件性,并分别记录两块残差。稠密求解可分解 m阶代数块及 n阶Schur块;成本还含多右端求解和矩阵乘法,直接组装的一轮为 O((n+m)3)量级,函数及导数评值另计。

把一个离散步完整复算 ​

对式(7)取 h=1/4,式(9)成为

5X=4+Z,Z2=X+14.

消元得到 20Z2−4Z−21=0。与 Z=1连续相接的根为

(11)Z=1+10610,X=41+10650.

另一根落在负分支,不是本初值的小步延续。将式(11)代回两组残差,它们都精确为零;这认证了后向Euler步,不代表已经得到原DAE在 t=1/4的精确轨道值。

大型模型还可能给出奇异质量矩阵 Ex′=Ax+f。此时正则常系数DAE会区分能够自由设定的微分初值和被外力及其导数强制确定的代数分量;单凭某一步矩阵可逆,仍不足以判断输入初值一致。

参考资料
  • Ernst Hairer,Solving Differential Equations on Manifolds,June2011,IV.2,印刷pp.30–32;V.2,p.42起:半显式指数1系统、隐函数约化及数值法与约化的关系。本文显式保留时间偏导 gt,并给出两块残差和具体分支计算。
  • Steffen Schulz,Four Lectures on Differential-Algebraic Equations,2003,Introduction及§2.2:隐式微分关系、导数缺失与指标。即时残差、一致状态和一致导数在本页分别定义。
关系图谱16 个相邻概念 · 2 类关系

拖动节点调整位置。

显示关系

显示:依赖

  1. 前置三跳
  2. 前置二跳
  3. 前置一跳
  4. 当前条目
  5. 后续一跳
  6. 后续二跳
  7. 后续三跳
文字版关系按与当前条目的最短距离分组
类型化关系