Constrained dynamics shake and rattle

For each parameter multi-index $\beta$ with $\mathrm{\mid }\beta \mathrm{\mid }\le k$, iteratively apply the chain rule [F2] to the identities ${\mathrm{\partial }}_h^{j}f(0,\alpha )=0$. Thus ${\mathrm{\partial }}_{\alpha }^{\beta }{\mathrm{\partial }}_h^{j}f(0,\alpha )=0$ for every $0\le j<r$, as needed below.

For each parameter multi-index $\beta$ with $\mathrm{\mid }\beta \mathrm{\mid }\le k$, put $k_{\beta }=k-\mathrm{\mid }\beta \mathrm{\mid }$, $N=r+k_{\beta }$, and $g_{\beta }(h)={\mathrm{\partial }}_{\alpha }^{\beta }f(h,\alpha )$. For every $s\le k_{\beta }$, put $m=N-s$. If $m\ge 2$, apply [F1] componentwise to $g_{\beta }^{(s)}$ through degree $m-1$; continuity of its highest derivative lets us replace the Lagrange coefficient at $\theta h$ by its value at zero plus $o(1)$. If $m=1$, subtract the linear Taylor polynomial: the derivative of the remainder is $o(1)$ by continuity, so [F3] gives a remainder $o(\mathrm{\mid }h\mathrm{\mid })$. Thus in every case $g_{\beta }^{(s)}(h)$ equals its Taylor polynomial through degree $m$ plus $o(\mathrm{\mid }h{\mathrm{\mid }}^m)$. By step 1.1 all coefficients of degrees below $r$ vanish. For $h\mathrm{\ne }0$, the quotient rule expresses the $\mathrm{\ell }$th derivative of $g_{\beta }\mathrm{/}{h}^r$ as a finite linear combination of $h^{-r-\mathrm{\ell }+s}{g}_{\beta }^{(s)}(h)$, $0\le s\le \mathrm{\ell }$. Substituting the expansions shows that the remainder terms are $o(\mathrm{\mid }h{\mathrm{\mid }}^{k_{\beta }-\mathrm{\ell }})$; the polynomial part has a continuous extension whose $\mathrm{\ell }$th derivative at zero is $\mathrm{\ell }!{\mathrm{\partial }}_h^{r+\mathrm{\ell }}{\mathrm{\partial }}_{\alpha }^{\beta }f(0,\alpha )\mathrm{/}(r+\mathrm{\ell })!$. Continuity of the derivatives on compact parameter neighborhoods makes these limits locally uniform in $\alpha$.

Define the value at $h=0$ by the displayed limit formula. For each $\mathrm{\mid }\beta \mathrm{\mid }+\mathrm{\ell }\le k$, step 2.1 gives a locally uniform limit for the corresponding derivative as $h\to 0$; at $h=0$ these formulas are differentiable in $\alpha$ because $f$ is $C^{k+r}$, and their $h$-derivatives agree with the next displayed coefficient. Induction on total derivative order therefore shows that all mixed partial derivatives through order $k$ exist and are continuous across $h=0$. The extension is $C^k$ with the asserted formula.

At any $q\in Q_{\mathrm{r}\mathrm{e}\mathrm{g}}$, full row rank of $J(q)$ gives a permutation of coordinates $q=(x,y)$ for which the $m\times m$ block $D_{y}c(q)$ is invertible. Apply [F2] to the smooth map $c$ for each finite $k\ge 1$. The resulting local graphs agree with the unique $C^1$ graph on overlaps, so that graph is $C^k$ locally for every $k$ and hence smooth. Thus locally $Q_{\mathrm{r}\mathrm{e}\mathrm{g}}$ is $y=\psi (x)$; $\varphi (x)=(x,\psi (x))$ parametrizes it, and $D\varphi$ has rank $d-m$ with image $\ker J$.

Define $G(x)=D\varphi (x{)}^{T}MD\varphi (x)$. For every nonzero $\xi$, $D\varphi \xi \mathrm{\ne }0$ and ${\xi }^{T}G(x)\xi =(D\varphi \xi {)}^{T}M(D\varphi \xi )>0$, so $G$ is invertible. Its entries are smooth, and $G^{-1}=\mathrm{adj}(G)\mathrm{/}\det G$ is smooth because $\det G>0$. For a cotangent coordinate $\pi \in {\mathbb{R}}^{d-m}$, put $p=MD\varphi G^{-1}\pi$. Then $J{M}^{-1}p=JD\varphi G^{-1}\pi =0$ and $D{\varphi }^{T}p=G{G}^{-1}\pi =\pi$. Conversely, if $J{M}^{-1}p=0$, then $M^{-1}p=D\varphi \xi$ for a unique $\xi$, and $D{\varphi }^{T}p=G\xi$; hence the displayed formula gives the unique ambient $p$ representing $\pi$. Both this map and its inverse $(q,p)\mapsto (x,D\varphi (x{)}^{T}p)$ are smooth, so $(x,\pi )\mapsto (\varphi (x),p)$ is a local diffeomorphism from $T^{\ast }{Q}_{\mathrm{r}\mathrm{e}\mathrm{g}}$ onto the part of $\mathcal{M}$ over this chart.

Let ${\theta }_0=p^{T}dq$ be the ambient canonical one-form. Under the chart just constructed, $dq=D\varphi dx$ and $D{\varphi }^{T}p=\pi$, so its pullback is ${\pi }^{T}dx$. Therefore the pullback of ${\omega }_0=-d{\theta }_0$ is ${\sum }_{j=1}^{d-m}d{x}_j\wedge d{\pi }_j$, which is closed and nondegenerate.

Every point of $\mathcal{M}$ has a chart of this form; the local calculation proves the restricted two-form is closed and nondegenerate everywhere. By [F3] it is symplectic, as claimed.

At a point of the compact regular chart $K\subset Q_{\mathrm{r}\mathrm{e}\mathrm{g}}$, permute coordinates as $q=(x,y)$ so $D_{y}c$ is an invertible $m\times m$ block. Apply [F2] for each finite regularity order; uniqueness makes these graphs agree locally, so the graph $y=\psi (x)$ is smooth. Write $\varphi (x)=(x,\psi (x))$ and $T(x)=D\varphi (x)$. Since $J(\varphi (x))T(x)=0$ and both kernels have dimension $d-m$, $\mathrm{im}T=\ker J$.

After shrinking to a convex coordinate neighborhood, $B(x)=D_{y}c(x,\psi (x))$ remains invertible and $\mathrm{\parallel }B(x{)}^{-1}(D_{y}c(x,y)-B(x))\mathrm{\parallel }\le \frac{1}{2}$. Set $R_x(y)=c(x,y)-B(x)(y-\psi (x))$. Then $R_x(\psi (x))=0$ and the derivative of $B(x{)}^{-1}{R}_x$ has norm at most $\frac{1}{2}$. The mean-value inequality [F4] gives $\mathrm{\parallel }B(x{)}^{-1}{R}_x(y)\mathrm{\parallel }\le \frac{1}{2}\mathrm{\parallel }y-\psi (x)\mathrm{\parallel }$. Since $B^{-1}c=(y-\psi )+B^{-1}{R}_x(y)$, it follows that $\mathrm{\parallel }q-\varphi (x)\mathrm{\parallel }=\mathrm{\parallel }y-\psi (x)\mathrm{\parallel }\le 2\mathrm{\parallel }B(x{)}^{-1}\mathrm{\parallel }\mathrm{\parallel }c(q)\mathrm{\parallel }$. Compactness makes the constant uniform on a finite chart cover.

For $\bar{q}=\varphi (x)$, put $G=T^{T}MT$ and define $\bar{p}=MT{G}^{-1}{T}^{T}p$. Positive definiteness of $M$ and injectivity of $T$ make $G$ invertible; its smooth entries and positive determinant give a smooth inverse by $G^{-1}=\mathrm{adj}(G)\mathrm{/}\det G$. Then $J(\bar{q})M^{-1}\bar{p}=0$ and $T^T(p-\bar{p})=0$. Since $\ker T^T$ and $\mathrm{im}J(\bar{q}{)}^T$ both have dimension $m$ and the latter is contained in the former, $p-\bar{p}=J(\bar{q}{)}^T\eta$ for a unique $\eta$.

The identity $J(\bar{q})M^{-1}(p-\bar{p})=A(\bar{q})\eta$ and [F1] imply $\mathrm{\parallel }p-\bar{p}\mathrm{\parallel }\le C\mathrm{\parallel }A(\bar{q}{)}^{-1}\mathrm{\parallel }(\mathrm{\parallel }J(q)M^{-1}p\mathrm{\parallel }+\mathrm{\parallel }J(\bar{q})-J(q)\mathrm{\parallel }\mathrm{\parallel }{M}^{-1}p\mathrm{\parallel })$. The momentum set is bounded and $J$ is locally Lipschitz, so step 2.1 bounds the right side by $C(\mathrm{\parallel }J(q)M^{-1}p\mathrm{\parallel }+\mathrm{\parallel }c(q)\mathrm{\parallel })$. Thus $\mathrm{\Pi }(q,p)=(\bar{q},\bar{p})$ has the asserted distance bound and lies in $\mathcal{M}$. If $(q,p)\in \mathcal{M}$, then $\bar{q}=q$ and the formula for $\bar{p}$ returns $p$, so $\mathrm{\Pi }$ is the identity there. By [F3] the coordinate formula is smooth, and it is the identity on the consistent phase set in this neighborhood. Shrink the finitely many neighborhoods so their closures lie in an open regular set; continuity keeps $A^{-1}$ and the required derivatives uniformly bounded there. A finite subcover gives finitely many such local retractions and a uniform constant.

Since $J(q)$ has full row rank and $M^{-1}$ is SPD, for every nonzero $\eta \in {\mathbb{R}}^m$ we have ${\eta }^{T}A(q)\eta =(J(q{)}^T\eta {)}^T{M}^{-1}(J(q{)}^T\eta )>0$. Thus $A(q)$ is invertible; its entries are smooth and $A^{-1}=\mathrm{adj}(A)\mathrm{/}\det A$, so its inverse and the displayed multiplier are smooth on a neighborhood of the chart.

Extend the formula in [F1] to an open neighborhood of $q_0$ on which $A$ remains invertible and define the ambient vector field $(q,p)\mapsto (M^{-1}p,-\mathrm{\nabla }V(q)-J(q{)}^T\lambda (q,p))$. It is smooth, hence locally Lipschitz, so [F3] gives a unique local solution through $(q_0,p_0)$. Shrink its time interval so $q(t)$ stays in that neighborhood.

Along this solution, the chain rule gives $\frac{d}{dt}c(q)=J(q)M^{-1}p$ and $\frac{d}{dt}(J(q)M^{-1}p)=D^{2}c(q)[v,v]-J(q)M^{-1}\mathrm{\nabla }V(q)-A(q)\lambda (q,p)=0$, where $v=M^{-1}p$. Since the second constraint starts at zero, it remains zero; the first derivative is therefore zero and $c(q)$ also remains zero. The chosen neighborhood is regular, and [F2] identifies its consistent constraint set as a smooth phase manifold, so the trajectory stays in $\mathcal{M}$.

The local uniqueness from [F3] applies to every solution of the same ambient extension, so in particular the constrained solution is unique.

Put $v=M^{-1}p$ and $w(\lambda )=-M^{-1}(\frac{1}{2}\mathrm{\nabla }V(q)+J(q{)}^T\lambda )$. Then $q_1=q+hv+h^{2}w(\lambda )$. Since $c(q)=0$ and $J(q)v=0$, [F3] extends $h^{-2}c(q_1)$ smoothly, with value $F_q(0,q,p,\lambda )=\frac{1}{2}{D}^{2}c(q)[v,v]-\frac{1}{2}J(q)M^{-1}\mathrm{\nabla }V(q)-A(q)\lambda$. Its $\lambda$ derivative is $-A(q)$, which is invertible. Since $A$ is smooth and invertible, $A^{-1}=\mathrm{adj}(A)\mathrm{/}\det A$ is smooth locally, so the center ${\lambda }_0$ varies smoothly with the phase input. Thus [F4] gives the unique local position multiplier branch at that center.

Let $G(h,q,p,\lambda ,\mu )=J(q_1)M^{-1}{p}_1$. Consistency gives $G(0,q,p,\lambda ,\mu )=0$, and [F3] extends $h^{-1}G$ smoothly. By [F5], its value at zero is $F_p(0,q,p,\lambda ,\mu )=D^{2}c(q)[v,v]-J(q)M^{-1}\mathrm{\nabla }V(q)-A(q)(\lambda +\mu )$. Its $\mu$ derivative is $-A(q)$. At $\lambda ={\lambda }_0$, the unique zero is $\mu ={\lambda }_0$, so [F4] gives a unique local momentum branch.

The branch equations are smooth in local phase coordinates, and $A^{-1}$ is uniformly bounded on the compact configuration projection of $C$. Since $J$ has full row rank on that compact projection, continuity gives an ambient neighborhood where it remains full row rank. Around each point of $C$, shrink the phase and multiplier neighborhoods so both roots remain in common IFT uniqueness neighborhoods on overlaps and exact $q_1$ remains in this regular neighborhood. These neighborhoods cover $C$; compactness gives a finite subcover and the minimum of their positive step radii is a common $h_0(C)$. The smooth branch derivatives are bounded on smaller compact neighborhoods covering $C$ for $\mathrm{\mid }h\mathrm{\mid }\le h_0(C)$, and continuity of the multiplier Jacobians preserves uniform invertibility after a further shrink.

The exact first constraint equation gives $c(q_1)=0$, the exact final projection equation gives $J(q_1)M^{-1}{p}_1=0$, and [F3] keeps $q_1$ in $Q_{\mathrm{r}\mathrm{e}\mathrm{g}}$. By [F1], this is exactly $(q_1,p_1)\in \mathcal{M}$.

Starting from the forward equations, use the reverse step size $-h$ and choose its first multiplier to be the forward $\mu$. Its first half-momentum is then $p_1+\frac{h}{2}\mathrm{\nabla }V(q_1)+hJ(q_1{)}^T\mu =p_{1\mathrm{/}2}$, so its next position is $q_1-h{M}^{-1}{p}_{1\mathrm{/}2}=q$.

Choose the reverse step's final multiplier to be the forward $\lambda$. Its final momentum is then $p_{1\mathrm{/}2}+\frac{h}{2}\mathrm{\nabla }V(q)+hJ(q{)}^T\lambda =p$. The reverse constraint equations hold because the endpoints are the forward consistent states. As $h\to 0$, the swapped multipliers and the reverse input converge to the same regular-chart branch centers; after shrinking $h_0$, both multipliers lie in the reverse IFT neighborhoods. By uniqueness of the selected branches [F3], this reverse solve is the defined $-h$ step. Thus ${\mathrm{\Psi }}_{-h}{\mathrm{\Psi }}_h=\mathrm{i}\mathrm{d}$; interchanging the endpoints and the multipliers proves ${\mathrm{\Psi }}_h{\mathrm{\Psi }}_{-h}=\mathrm{i}\mathrm{d}$.

If $h=0$, then ${\mathrm{\Psi }}_0$ is the identity and preserves ${\omega }_{\mathcal{M}}$. Assume $h\mathrm{\ne }0$. On a local configuration chart define $L_h(q,q_1)=\frac{1}{2h}(q_1-q{)}^{T}M(q_1-q)-\frac{h}{2}(V(q)+V(q_1))$. Differentiating gives $-D_1{L}_h=p_{1\mathrm{/}2}+\frac{h}{2}\mathrm{\nabla }V(q)$ and $D_2{L}_h=p_{1\mathrm{/}2}-\frac{h}{2}\mathrm{\nabla }V(q_1)$. By [F3], the first expression differs from $p$ by the normal vector $-hJ(q{)}^T\lambda$, and the second differs from $p_1$ by $hJ(q_1{)}^T\mu$. Each normal vector vanishes on tangent variations of $Q$, so the restricted cotangent covectors satisfy $p{\mathrm{\mid }}_{T_{q}Q}=-D_1{L}_h{\mathrm{\mid }}_{T_{q}Q}$ and $p_1{\mathrm{\mid }}_{T_{q_1}Q}=D_2{L}_h{\mathrm{\mid }}_{T_{q_1}Q}$.

Let $\theta$ be the canonical one-form on $T^{\ast }Q$, ${\theta }_{(q,\pi )}(\delta q,\delta \pi )=\pi (\delta q)$. The identities in step 1.1 give ${\mathrm{\Psi }}_h^{\ast }\theta -\theta =d(L_h(q,q_1(q,p)))$ on the input chart, because both sides applied to an input variation equal $D_1{L}_h\delta q+D_2{L}_h\delta q_1$.

Taking the exterior derivative of step 2.1 yields ${\mathrm{\Psi }}_h^{\ast }d\theta =d\theta$. By [F1], ${\omega }_{\mathcal{M}}$ is the restriction of the ambient canonical form; [F2] identifies it as $-d\theta$ in cotangent coordinates. Hence this is ${\mathrm{\Psi }}_h^{\ast }{\omega }_{\mathcal{M}}={\omega }_{\mathcal{M}}$. The branch map is smooth by [F4] and locally invertible by [F6]. By [F2] and [F5], the preserved form is the symplectic form on $\mathcal{M}$; this proves the stated symplecticity.

If $h=0$, [F1] returns exact-identity immediately. For nonzero $h$, consider either chord call with budget $N$ and refresh period $1$. Its initial iterate lies in its open domain. At index $k$, its residual evaluation and norm comparison terminate by the admitted primitive contract. The tolerance and budget branches return the current valid iterate and its actual residual norm. If neither branch returns, then $k<N$; one terminating Jacobian evaluation and factorization attempt either return singular-jacobian with this same data or supply a retained factor. Its terminating exact solve produces the correction. The returning domain test either returns invalid-trial with the old valid iterate and evaluated residual or admits the new iterate and increments the index. Thus induction preserves domain membership and validity of all returned residual data.

Every continuing pass strictly increments $k$, and an entered pass $k=N$ returns before a factorization or correction. Consequently there are at most $N+1$ residual checks and at most $N$ correction attempts, Jacobian assemblies, factorization attempts, and retained-factor solves. Each pass consists of finitely many returning operations by step 1.1, so the call terminates, including when $N=0$. These are counts of loop operations, not a uniform bound on the internal cost of a supplied evaluation or backend.

Apply step 2.1 first with $N=N_q$. A non-root position status is returned immediately with its last valid multiplier and actual residual, without a momentum call. On position success the finite assignments in [F1] form the momentum residual. The all-$\lambda \in U_q$ admissibility condition ensures that even this approximate returned position multiplier supplies a $C^1$ momentum residual on the whole open $U_p$ with returning evaluations. Apply step 2.1 with $N=N_p$ to the second call. Any failure returns the last valid momentum multiplier and actual tangency residual together with the successful position-stage values. If it succeeds, the remaining assignments and output norms terminate and return the prescribed state, both multipliers, residual norms, and success status. There are at most two calls, at most $N_q+N_p+2$ loop residual checks, and at most $N_q+N_p$ correction/factorization/solve attempts in total. No branch or count asserts success for every finite budget, an exact zero from a positive tolerance, or a forward-error certificate.