Skip to content

交付一份DAE相容、指标与漂移证书 ​

微分代数系统的相容与漂移路线要求最终交付三份可以独立复算的记录。它们分别回答:这个初值能否产生轨道,哪些状态其实由数据导数强制决定,以及反馈后的约束误差能保证到什么程度。

若还不熟悉行列式,先复习它对可逆性的判据;任务一用它检查代数块,任务二用它判定矩阵束是否正则。任务三还使用对称正定矩阵作为质量矩阵,并据此认证约束乘子的唯一性。

下载公开复算程序及结果。精确矩阵、有理不等式与多项式部分用Fraction核验,高精度计算用于复核导数、卷积和显示小数;一般结论仍由正文证明负责。

任务一:两个耦合代数变量不能分别随意初始化 ​

考虑

x′=−x+z1,z1+z2−x−t=0,z1z2−1=0,(x,z1,z2)(0)=(5/2,2,1/2).

交出一致导数,再取后向Euler步长 h=1/2,给出与初值连续相接的完整三元组。验收须同时包括两项代数残差、微分残差和分支判断,不能只报一个非线性求解器的成功标志。

一致导数 ​

代数Jacobian及行列式为

gz=(11z2z1),det⁡gz=z1−z2=3/2.

所以可逆代数块定理适用。先算 x0′=−1/2,再微分两条约束:

z1′+z2′=x0′+1=1/2,12z1′+2z2′=0.

解得

(x0′,z1′,z2′)=(−1/2,2/3,−1/6).

显式时间的导数1不可丢掉;若误把约束当成自主系统,会算出另一组不一致导数。

隐式一步与根的连续选择 ​

后向Euler方程化为

3X=5+Z1,Z1+Z2=X+1/2,Z1Z2=1.

消去 X,Z2,得到 4Z12−13Z1+6=0。所需解是

(1)Z1=13+738,Z2=13−7312,X=53+7324.

其中乘积由 (169−73)/96=1直接检查,其余两项也可逐项消去。

要认证分支,可让步长从0增加到1/2。一般消元式是

Z12−(5/2+h+h2)Z1+(1+h)=0.

其判别式在这段区间严格为正:在0处等于9/4,导数为 2(5/2+h+h2)(1+2h)−4≥1。大根从2连续出发,小根从1/2出发;式(1)取大根,对应原始 z1>z2分支。

实际三方程Newton矩阵为

Jh=(1+h−h0−1110Z2Z1).

在式(1)处,其行列式是 Z1−(3/2)Z2=73/4>0。这是该根局部非退化的证书;它没有保证任意Newton初猜都到达大根,也没有把离散步等同于原轨道值。

任务二:混合坐标后找回自由初值和三阶链 ​

给定 Ex′=Ax+f(t),其中

E=(1−11−1−11001−101−110−1),A=(−11−111000−101−110−12),f(t)=(t−ttt3−t).

不能从 E哪一行有零元素来猜微分变量。使用下面的左右变换证书:

P=(1000110001100011),Q=(1100011000110001).

两矩阵行列式均为1。精确相乘应得到

(2)PEQ=diag(1,N),PAQ=diag(−1,I3),N=(010001000),Pf=(t,0,0,t3)T.

全部一致状态 ​

写 x=Q(y,z1,z2,z3)T。由正则常系数分解,

y′=−y+t,z2′=z1,z3′=z2,0=z3+t3.

因此 z3=−t3,z2=−3t2,z1=−6t,只有 y(0)可以自由选择。取 y(0)=2,得到

(3)y=t−1+3e−t,x(t)=(−5t−1+3e−t−6t−3t2−3t2−t3−t3).

在 t=0,一致状态必须形如 (c,0,0,0);所选状态为 (2,0,0,0),一致导数为 (−8,−6,0,0)。直接代回 Ex′(0)=Ax(0)+f(0)只能验证一层必要条件,完整相容性来自式(2)的整条链。

束行列式为 −(λ+1),并非恒零。又有 N3=0,N2≠0,所以在保留这四个状态分量的表示下,微分指标为3。两个代数微分已经决定 z1,第三个才决定 z1′。

微小外力如何放大成很大的状态 ​

把原坐标的外力改成

fω=f+(0,0,0,ω−1sin⁡(ωt))T,ω>0,

并保留同一个自由初值 y(0)=2。由于 P的最后一列仍为 (0,0,0,1)T,变动只进入链末端。各代数变化准确为

Δz3=−ω−1sin⁡(ωt),Δz2=−cos⁡(ωt),Δz1=ωsin⁡(ωt).

所以 ‖fω−f‖∞≤1/ω,但原坐标第一分量有 Δx1=ωsin⁡(ωt)。当 ω≥π/2,在时间 t=π/(2ω)∈[0,1]已经达到幅度 ω。

扰动后的强制初值也随之改变:Δx(0)=(0,−1,−1,0)。这里没有假装两份不同强制外力仍能共享全部原始初值;共享的只有真正自由的 y(0)。证书展示的是数据微分损失,控制数据本身的幅度不够,还要控制系统所调用的导数。

任务三:非单位质量椭圆上的反馈与容差 ​

取

M=diag(2,3),F=(2,−1)T,g(q)=q12+2q22−12.

采用动力符号 Mv′=F−GTλ,其中 G=(q1,2q2)。这是双侧等式约束,乘子不另受单侧接触力的符号条件限制。

先认证原约束状态 ​

对

q=(1/3,2/3),v=(−4/5,1/5),

有 g=0,Gv=0。又

H=v12+2v22=18/25,W=GM−1GT=35/54,GM−1F=−1/9.

因此原约束乘子与加速度为

λ=822875,v′=(738875,−657875)T.

逐项核验 Mv′+GTλ=F和 Gv′+H=0,即得到动力与加速度约束的双证书。它使用真实质量度量,不能把 W替换成 GGT。

从不一致位置开始,反馈究竟承诺什么 ​

现在把位置改为 q∗=(2/5,4/5),速度仍取上面的 v。得到

e0=g(q∗)=11/50,w0=G(q∗)v=0,W∗=14/15.

此状态不属于原DAE的允许初值。取临界阻尼反馈 α=β=3,则

(4)λ∗=114,v∗′=(920,−95)T.

核验时使用 G∗=(2/5,8/5)、H=18/25:动力行仍精确成立,而

G∗v∗′+H=−99/50=−9e0.

如果不加反馈,只令加速度约束残差为零,则 e″=0,这个初始偏差会恒等于11/50。精确稳定化轨道则在其有效存在区间满足

(5)e(t)=1150(1+3t)e−3t.

反馈逐渐修正偏差,没有在 t=0把它投回椭圆。

一份带持续缺陷的有限时间证书 ​

设已获得定义在 [0,2]上的 C2候选位置曲线,取 v=q′,其初始残差为上述 e0,w0,并已独立认证

|η(t)|=|e″+6e′+9e|≤9100000(0≤t≤2).

临界阻尼核界给

(6)|e(2)|≤7750e−6+1100000(1−7e−6)=1100000+153993100000e−6<1250.

最后的不等式可以完全用有理数核验:

e6>∑k=0136kk!=1005917325025>400,

再用 e−6<1/400代入式(6),得到上界 154393/40000000<1/250。这里既没有依赖显示小数,也没有假定持续缺陷为零。

该证书的输入是整个时间窗上的曲线与缺陷界。只有若干离散样本的残差小,不能直接升级为这个连续合同;还须给出插值及区间内缺陷控制。

连续衰减与显式失稳可以同时出现 ​

对残差系统 (e,w)′=(w,−9e−6w),令 u0=(e0,0),用显式Euler步长 h=2/3。更新矩阵的双特征值为−1,且不是对角块,直接得到

en=e0(−1)n−1(2n−1),n≥1.

第四步已经给 e4=−77/50,而连续解式(5)一直为正并衰减。若选 h=1/10,则 hω=3/10<2,该线性测试的每个初值都渐近衰减。

这项离散测试没有额外声称完整椭圆积分器稳定。对物理位置本身作Euler步,准确展开是

g(q+hv)=g(q)+hG(q)v+h22(v12+2v22).

即使从前面的相容状态出发,位置步仍增加 9h2/25。精确连续反馈、稳定步长和逐步精确保持约束,必须分别验收。

最终交付 ​

交付原方程和所有变量的状态角色、允许初值及其导数、使用的可逆块、矩阵束变换、全部隐藏约束和剩余自由参数。数值部分另列步长、代数分支、各类残差、持续缺陷的时间窗,以及究竟认证了离散步还是连续曲线。

上述三份证书分别改变了代数耦合、坐标表示和质量/约束几何。若更换约束或求解器,应重新推导这些接口;单独沿用某个乘子数值或某个步长阈值,不能保留原来的保证。