Baumgarte约束稳定化 | Codeboy's Blog

把位置约束求两次导数，可以算出反作用力，却也会丢掉两个积分常数。求解器若只推进加速度层，微小的初始位置和速度误差便可能变成长期漂移。Baumgarte方法保留已经偏离了多少、偏离得多快这两项信息，再把它们反馈给加速度。

形式陈述

在独立位置约束附近定义反馈

取常对称正定矩阵 $M\in\mathbb R^{d\times d}$ 、光滑外力 $F(t,q,v)$ 及光滑约束 $g:U\to\mathbb R^m$ ，其中 $U\subset\mathbb R^d$ 开、 $1\le m\le d$ 。假设 $G(q)=Dg(q)$ 在所讨论邻域满行秩。原位置约束DAE是

q'=v,\qquad Mv'=F-G^T\lambda,\qquad g(q)=0. \tag{1}

记

H(q,v)=D^2g(q)[v,v],\qquad W(q)=GM^{-1}G^T.

$H$ 是一个 $m$ 维向量，其第 $i$ 项为 $v^TD^2g_i(q)v$ 。取常数 $\alpha>0,\beta>0$ ，将乘子定义为

\lambda=W^{-1}\bigl(GM^{-1}F+H+2\alpha Gv+\beta^2g\bigr). \tag{2}

式(2)与前两条运动方程组成Baumgarte稳定化系统。它在约束面附近有定义，也允许从轻微偏离 $g=0$ 的状态开始。这里并没有继续另加一个不相容的代数要求 $g(q)=0$ 。

连续残差的精确合同

沿稳定化系统的任意经典解，令 $e(t)=g(q(t))$ ，则

e''+2\alpha e'+\beta^2e=0. \tag{3}

因此，在解存在且约束保持独立的时间段内，残差按这个常系数方程演化。如果 $e(0)=e'(0)=0$ ，则 $e$ 恒为零，稳定化轨道也是原约束系统的轨道。

若初值不相容，反馈使残差趋零的结论需要轨道一直存在于有效区域；它不会将非零初始残差立即改成零，也不能使该初值成为原DAE的一致初值。

直觉

为什么这个乘子能控制漂移

由于 $G$ 满行秩，任意非零 $a\in\mathbb R^m$ 都有 $G^Ta\ne0$ 。由正定性，

a^TWa=(G^Ta)^TM^{-1}(G^Ta)>0,

所以 $W$ 可逆，式(2)确实唯一确定乘子。链式法则给

e'=Gv,\qquad e''=Gv'+H=GM^{-1}F-W\lambda+H.

代入式(2)，所有原动力项抵消，剩下式(3)。这说明两个反馈项的符号为什么必须是当前约定中的正号；若动力方程改用 $+G^T\lambda$ ，乘子的符号也要一起改变。

不加反馈时只规定 $e''=0$ ，会留下 $e(t)=e(0)+te'(0)$ 。加入 $2\alpha e'$ 抑制残差变化，加入 $\beta^2e$ 消除静态偏差，两者负责的错误不同。

衰减来自哪里

对每个残差分量，特征根是

-\alpha\pm\sqrt{\alpha^2-\beta^2}.

当 $\alpha,\beta>0$ 时，两根实部都严格为负：实根情形有 $\sqrt{\alpha^2-\beta^2}<\alpha$ ，共轭根情形实部就是 $-\alpha$ 。由常系数线性ODE的解形式，残差与残差速度均衰减；重根时还会乘一次多项式因子。

$\alpha=0,\beta>0$ 只给振荡， $\beta=0,\alpha>0$ 允许非零常数残差。因此“至少放一个很大的反馈系数”并不是两类漂移都消失的条件。

一致初值为何保留原轨道

如果 $g(q_0)=0$ 且 $G(q_0)v_0=0$ ，式(3)的零初值唯一解就是零。此时反馈项在整个轨道上消失，式(2)退回原系统由二阶隐藏约束确定的乘子。反过来，原系统的每条光滑轨道也满足稳定化系统。

这个等价只针对相容初值和同一有效邻域。约束梯度失秩、解逃离邻域或发生有限时间爆破时，线性残差方程不能单独保证物理轨道继续存在。

例子与边界

临界阻尼给出直接可用的误差界

取 $\alpha=\beta=\omega>0$ 。更一般地，对一个实际 $C^2$ 位置曲线 $q$ ，令 $v=q'$ ，并定义它的约束方程缺陷为

\eta=e''+2\omega e'+\omega^2e. \tag{4}

在 $0\le t\le T$ 上，常数变易得到

e(t)=\bigl[e_0+(w_0+\omega e_0)t\bigr]e^{-\omega t} +\int_0^t(t-s)e^{-\omega(t-s)}\eta(s)\,ds, \quad w_0=e'(0). \tag{5}

它可以直接求两次导数并代回式(4)检查，不需要把计算出的轨道先当成精确解。

对任意向量范数，若 $\sup_{[0,T]}\|\eta\|\le\delta$ ，则

\begin{aligned} \|e(t)\|\le{}& \bigl(\|e_0\|+t\|w_0+\omega e_0\|\bigr)e^{-\omega t}\\ &+\frac{\delta}{\omega^2} \bigl[1-(1+\omega t)e^{-\omega t}\bigr]. \end{aligned} \tag{6}

最后一项来自非负核的准确积分。若该缺陷界在所有正时间成立，相同公式给全时间界；若只在 $[0,T]$ 核验，就只在该窗使用。再用 $te^{-\omega t}\le1/(e\omega)$ 可得较粗的统一界

\|e(t)\|\le\|e_0\|+ \frac{\|w_0+\omega e_0\|}{e\omega}+ \frac{\delta}{\omega^2}. \tag{7}

例如 $e_0=1/100$ 、 $w_0=0$ 、 $\omega=2$ 、缺陷恒为 $\eta=1/250$ 时，

e(t)=\frac1{1000}+\frac9{1000}(1+2t)e^{-2t}. \tag{8}

极限为 $1/1000$ ，准确达到持续缺陷的静态尺度 $\delta/\omega^2$ 。稳定并不等于误差最后必为零。

增大反馈会给显式步长增加限制

先只看线性残差测试系统

\frac{d}{dt}\binom e w=B\binom e w,\qquad B=\begin{pmatrix}0&1\\-\omega^2&-2\omega\end{pmatrix}. \tag{9}

对它使用显式Euler，更新矩阵是 $R=I+hB$ ，其双特征值为 $r=1-h\omega$ 。由绝对稳定域，严格衰减要求

0<h\omega<2. \tag{10}

重根还要检查Jordan因子。写 $C=B+\omega I$ ，有 $C^2=0$ 且 $C\ne0$ ，于是对 $r\ne0$ 、整数 $n\ge1$ ，

R^n=r^nI+nh\,r^{n-1}C. \tag{11}

当 $|r|<1$ ，两项都趋零； $r=0$ 时直接有 $R^2=0$ 。边界 $h\omega=2$ 给 $r=-1$ ，对满足 $Cu_0\ne0$ 的初值出现线性增长；超过边界时可沿特征方向指数增长。连续方程越来越快的衰减，可能在固定显式步长下变成数值不稳定。

这个二阶测试只诊断反馈引入的快尺度，不独自证明完整非线性机械离散法稳定。特别是，对 $q,v$ 作一步Euler后再算非线性 $g(q)$ ，所得残差一般不等于直接对式(9)做Euler。

左图的零缺陷连续解展示反馈对初始漂移的修正；右图固定残差初态，用步数比较三个归一化步长，边界重根仍会增长。两侧计算对象不同，右图没有将非线性位置残差直接当成线性Euler递推。

精确连续约束也不等于逐步精确约束

取 $g(q)=(\|q\|^2-1)/2$ ，一步位置更新 $q_+=q+hv$ 给

g(q_+)=g(q)+hG(q)v+\frac{h^2}{2}\|v\|^2. \tag{12}

即使当前 $g(q)=G(q)v=0$ ，只要 $v\ne0$ ，该显式位置步仍会离开圆面。这与式(3)的连续结论没有矛盾：离散点并不是该微分方程的精确轨道。

RATTLE在步内另求位置与切向动量约束。Baumgarte反馈本身没有这两次投影，也不自动获得保辛、保能量或逐步零约束误差。

推论与应用

将反馈装进一个可核验的乘子求解

一次加速度评值可直接解

\begin{pmatrix}M&G^T\\G&0\end{pmatrix} \binom a\lambda= \binom F{-H-2\alpha Gv-\beta^2g}. \tag{13}

消去 $a=M^{-1}(F-G^T\lambda)$ ，就恢复式(2)。验证程序应同时记录这两行残差，以及位置残差 $g$ 和速度残差 $Gv$ 。

对满足 $q'=v$ 的光滑候选曲线，若第二行残差定义为

r_2=Gv'+H+2\alpha Gv+\beta^2g,

那么 $r_2$ 就是式(4)的 $\eta$ ，可直接进入式(6)。仅知道动力第一行残差小，并不能控制 $\eta$ 。若插值曲线还不满足 $q'=v$ ，链式法则会多出运动学缺陷项，必须先把它们算入；不能把任意离散求解器的一项残差直接冒充连续误差证书。

稠密常质量矩阵可以预分解。一次评值形成 $W$ 、应用质量矩阵逆并解 $m$ 维正定系统；朴素稠密算术成本为 $O(d^2m+dm^2+m^3)$ ，另加 $F,G,H$ 的评值。这里应用预分解的逆指三角求解，不要求显式储存 $M^{-1}$ 。约束接近相关时， $W$ 会病态，反馈系数调大不能恢复丢失的独立约束。

从连续容差走到离散设置

先声明约束和速度残差采用的单位及尺度，再按允许的衰减时间选 $\alpha,\beta$ 。临界阻尼下，式(6)把初始偏差、持续缺陷和时间窗分开；随后检查积分器能否解析 $1/\omega$ 这一快尺度。

最终验收应交出原初值是否一致、反馈方程是否被正确实现、离散步长是否稳定，以及实际位置/速度约束误差各有多大。四个判断分别需要证据，不能由一次成功的乘子线性求解代替。

参考资料

J. Baumgarte，Stabilization of constraints and integrals of motion in dynamical systems，Computer Methods in Applied Mechanics and Engineering 1(1)，1972，pp.1–16：方法的原始论文；本页符号与误差界在正文独立推导。
RBDL官方文档，Constraints / Baumgarte stabilization，反馈项 $-2\alpha\dot\phi-\beta^2\phi$ 与加速度乘子系统；Constraint::setBaumgarteTimeConstant给出时间常数约定。本文固定自己的乘子符号并独立推导式(2)–(13)。
Ernst Hairer，Solving Differential Equations on Manifolds，June2011，IV.4，印刷pp.34–36：原位置、速度与加速度约束的关系。连续反馈界与离散重根边界在本文直接计算。

形式陈述 ​

在独立位置约束附近定义反馈 ​

连续残差的精确合同 ​

为什么这个乘子能控制漂移 ​

衰减来自哪里 ​

一致初值为何保留原轨道 ​

临界阻尼给出直接可用的误差界 ​

增大反馈会给显式步长增加限制 ​

精确连续约束也不等于逐步精确约束 ​

将反馈装进一个可核验的乘子求解 ​

从连续容差走到离散设置 ​

反向链接