Scaled RATTLE residuals bound the finite solve trajectory error

For $F_q(\lambda )=h^{-2}c(q_1(\lambda ))$, [F2] gives $D_{\lambda }{F}_q=-A(q)+o(1)$ uniformly on the branch ball. Shrink it so $\mathrm{\parallel }A(q{)}^{-1}(D_{\lambda }{F}_q+A(q))\mathrm{\parallel }\le \frac{1}{2}$. The derivative of $G(\lambda )=A(q{)}^{-1}{F}_q(\lambda )+\lambda$ is then bounded by $1\mathrm{/}2$. Since $F_q(\lambda )=0$ at the exact root $\lambda$, [F5] along the segment to $\tilde{\lambda }$ gives $\mathrm{\parallel }G(\tilde{\lambda })-G(\lambda )\mathrm{\parallel }\le \frac{1}{2}\mathrm{\parallel }\tilde{\lambda }-\lambda \mathrm{\parallel }$. As $A^{-1}{F}_q(\tilde{\lambda })=G(\tilde{\lambda })-G(\lambda )-(\tilde{\lambda }-\lambda )$, it follows that $\mathrm{\parallel }\tilde{\lambda }-\lambda \mathrm{\parallel }\le 2\mathrm{\parallel }A(q{)}^{-1}\mathrm{\parallel }{\epsilon }_q$. This inverse estimate is needed because [F4] alone does not convert residual size to state error.

The position formula gives ${\tilde{q}}_1-q_1=-h^2{M}^{-1}J(q{)}^T(\tilde{\lambda }-\lambda )$ and ${\tilde{p}}_{1\mathrm{/}2}-p_{1\mathrm{/}2}=-hJ(q{)}^T(\tilde{\lambda }-\lambda )$. The uniform matrix bounds and step 1.1 yield $\mathrm{\parallel }{\tilde{q}}_1-q_1\mathrm{\parallel }\le C{h}^2{\epsilon }_q$ and $\mathrm{\parallel }{\tilde{p}}_{1\mathrm{/}2}-p_{1\mathrm{/}2}\mathrm{\parallel }\le C\mathrm{\mid }h\mathrm{\mid }{\epsilon }_q$.

For fixed $\tilde{\lambda }$, the final scaled tangency residual is affine in $\mu$ with derivative exactly $-A({\tilde{q}}_1)$. Uniform invertibility bounds the distance from $\tilde{\mu }$ to the exact momentum root for this $\tilde{\lambda }$ by $C{\epsilon }_p$. The smooth momentum branch in [F2] changes that exact root by at most $C{\epsilon }_q$ when $\tilde{\lambda }$ replaces $\lambda$, using step 1.1 and the uniform branch derivative bound.

In the final momentum formula, the initial half-momentum difference contributes $O(\mathrm{\mid }h\mathrm{\mid }{\epsilon }_q)$ by step 2.1; changing the position in the smooth force and constraint Jacobian contributes $O(\mathrm{\mid }h{\mathrm{\mid }}^3{\epsilon }_q)$; changing the final multiplier contributes $O(\mathrm{\mid }h\mathrm{\mid }({\epsilon }_q+{\epsilon }_p))$ by step 2.2. Their sum proves $\mathrm{\parallel }{\tilde{p}}_1-p_1\mathrm{\parallel }\le C\mathrm{\mid }h\mathrm{\mid }({\epsilon }_q+{\epsilon }_p)$ and hence the one-step estimates.

The residual definitions give $\mathrm{\parallel }c({\tilde{q}}_1)\mathrm{\parallel }=\mathrm{\mid }h{\mathrm{\mid }}^2{\epsilon }_q$ and $\mathrm{\parallel }J({\tilde{q}}_1)M^{-1}{\tilde{p}}_1\mathrm{\parallel }=\mathrm{\mid }h\mathrm{\mid }{\epsilon }_p$. Whenever the output lies in a neighborhood $U_j$ from [F6], its local retraction changes it by at most $C(\mathrm{\mid }h{\mathrm{\mid }}^2{\epsilon }_q+\mathrm{\mid }h\mathrm{\mid }{\epsilon }_p)\le C\mathrm{\mid }h\mathrm{\mid }({\epsilon }_q+{\epsilon }_p)$ for $\mathrm{\mid }h\mathrm{\mid }\le 1$. For the repeated-trajectory claim, choose a compact phase neighborhood $D$ of ${\mathcal{C}}_1$ with closure inside the regular chart. The compact branch bounds and ${\mathrm{\Psi }}_0=\mathrm{i}\mathrm{d}$ imply ${\mathrm{\Psi }}_h({\mathcal{C}}_1)\subset D$ for sufficiently small $\mathrm{\mid }h\mathrm{\mid }$. Apply [F6] to $D$ and choose its finite retraction cover with a common distance constant. The union of these open domains contains the compact set $D$, hence contains some uniform ambient neighborhood of it. For every restored input in ${\mathcal{C}}_1$, step 3.1 puts its raw output within $O(\mathrm{\mid }h\mathrm{\mid }({\epsilon }_q+{\epsilon }_p))$ of the exact step from that same input in $D$. Bounded scaled residuals and sufficiently small $\mathrm{\mid }h\mathrm{\mid }$ therefore place the raw output in this union. The same-input exact step need not belong to the reference tube ${\mathcal{C}}_0$. Thus the retraction is defined, and the restored update differs from the exact RATTLE update by at most $C\mathrm{\mid }h\mathrm{\mid }({\epsilon }_q+{\epsilon }_p)$.

Work in the phase coordinates of the Statement and choose a slightly larger compact convex coordinate neighborhood on which the branches and their first derivatives are bounded by [F2]. Since ${\mathrm{\Psi }}_0$ is the identity, $D_z{\mathrm{\Psi }}_h=I+O(h)$ uniformly there, so [F5] gives a Lipschitz factor $1+Lh$ for small positive $h$. The one-step estimate and the retraction bound, with equivalence of ambient and chart norms on the compact tube, bound each restored solver perturbation by $Kh({\epsilon }_q+{\epsilon }_p)$. Starting with the same initial state, the error therefore obeys $e_{n+1}\le (1+Lh)e_n+Kh({\epsilon }_q+{\epsilon }_p)$ until a first exit from ${\mathcal{C}}_1$. Iteration gives $e_n\le K_T\epsilon$ for fixed residual bound and $e_n\le K_T{h}^2$ when each residual is at most $C_0{h}^2$. Choose $\epsilon$ small enough, or $h$ small enough in the second case, that this bound is less than the positive coordinate distance from ${\mathcal{C}}_0$ to the boundary of ${\mathcal{C}}_1$; then no first exit occurs. Step 4.1, applied until this first-exit argument closes, places every raw output in a retraction neighborhood for small $h$. Finally, the unscaled residuals are exactly $h^2{\epsilon }_q$ and $h{\epsilon }_p$, hence $O(h^4)$ and $O(h^3)$ in the second case.