= Accumulated contact-mode phase drift
{title2=$\Delta\phi=T(\widetilde\omega-\omega)$}
<Accumulated contact-mode phase drift> is the phase difference that builds over many timesteps when a numerical contact mode oscillates at a slightly incorrect <frequency>. A small error per step can matter after many cycles even when the amplitude remains bounded. The fast oscillations arise naturally in stiff compliant-contact models, whose numerical integration difficulty is described in the primary multibody notes.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, Compliant Contact Models]{ref}
For a fixed undamped mode of angular <frequency> $\omega$, step $h$, and numerical phase advance $\theta(h\omega)$ per step, after $N$ steps the unwrapped phase error is
$$
\boxed{\Delta\phi_N=N[\theta(h\omega)-h\omega]
=T(\widetilde\omega-\omega),\qquad T=Nh.}
$$
Unwrapped phase counts complete cycles; reducing it modulo $2\pi$ would hide a missed or added cycle. This expression assumes the same fixed mode and a uniform step throughout the comparison.
For central differences, $\theta(z)=2\arcsin(z/2)$ on $0<z<2$. Expanding gives
$$
\Delta\phi_N=\frac{T\omega^3h^2}{24}+O(T\omega^5h^4).
$$
For the trapezoidal rule, $\theta(z)=2\arctan(z/2)$ instead gives
$$
\Delta\phi_N=-\frac{T\omega^3h^2}{12}+O(T\omega^5h^4).
$$
Both are second-order phase errors at fixed duration, but their signs and leading coefficients differ. Halving $h$ approximately quarters the accumulated phase error; doubling the simulated duration approximately doubles it when the mode remains unchanged.
For an original example, choose $\omega=600\,\mathrm{rad/s}$, $h=0.5\,\mathrm{ms}$ and $T=0.20\,\mathrm s$, so $N=400$ and $h\omega=0.3$. Central differences give $\widetilde\omega\approx602.273\,\mathrm{rad/s}$ and a phase lead about $0.4546\,\mathrm{rad}$, or $26.05^\circ$. The relative <frequency> error is only about $0.379\%$, yet the phase lead has become readily visible. The leading-order estimate is $0.4500\,\mathrm{rad}$.
An approximate phase budget $|\Delta\phi_N|\le\phi_{\max}$ for central differences suggests
$$
h\lesssim\sqrt{\frac{24\phi_{\max}}{T\omega^3}}.
$$
With the preceding mode and duration and $\phi_{\max}=0.020\,\mathrm{rad}$, it gives $h\lesssim105.4\,\mathrm{\mu s}$. Taking $100\,\mathrm{\mu s}$ gives an exact phase lead about $0.01801\,\mathrm{rad}$. This is far below the $3.333\,\mathrm{ms}$ oscillator stability boundary; <contact-mode timestep stability> and timing accuracy impose different requirements.
The associated timing shift is approximately $-\Delta\phi/\omega$ for a fixed phase feature of a nearly sinusoidal mode with small drift. It is not a general contact-event error bound: a changing amplitude, a grazing zero, changing tangent <stiffness> or an active-contact transition can invalidate that approximation. A variable-step method also needs its actual update analysis; summing uniform-step formulas without checking how successive maps compose is not automatically exact.
A proposed convergence test records a small compressed-contact oscillation over a stated duration, compares unwrapped numerical and reference phases, and repeats with half the step and twice the duration. Compare the <force> waveform and event times as well as the maximum penetration. Material damping changes the physical response, while numerical phase drift changes the timing of that modeled response; they should be tested separately.
Back to article page
