= Contact-mode numerical phase error
<Contact-mode numerical phase error> is the difference between the oscillation phase produced by a time integrator and the phase of the modeled contact dynamics. It can grow while the numerical solution remains bounded. <Contact-mode timestep stability> controls a different property; stable integration alone does not ensure accurate <force> timing or contact-event timing. The physical motivation is the fast oscillation introduced by a stiff compliant contact model.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, Compliant Contact Models, stiff-spring and numerical integration discussion]{ref}
For a frozen undamped contact mode $\ddot x+\omega^2x=0$, central differences give
$$
x_{n+1}-2x_n+x_{n-1}=-h^2\omega^2x_n.
$$
Put $x_n=e^{in\theta}$ and use $e^{i\theta}+e^{-i\theta}=2\cos\theta$. On the stable branch $0<h\omega<2$,
$$
\cos\theta=1-h^2\omega^2/2,\qquad
\boxed{\widetilde\omega=\frac{\theta}{h}=\frac2h\arcsin(h\omega/2).}
$$
The numerical <frequency> exceeds the modeled <frequency>. A <Taylor series> gives
$$
\widetilde\omega=\omega+\frac{\omega^3h^2}{24}+O(\omega^5h^4).
$$
For the trapezoidal rule, its amplification-factor phase instead gives
$$
\widetilde\omega=\frac2h\arctan(h\omega/2)
=\omega-\frac{\omega^3h^2}{12}+O(\omega^5h^4).
$$
That scheme preserves modal amplitude in this undamped linear setting while lagging in phase. These formulas isolate a single mode with fixed <frequency> and consistent initialization; general nonlinear contacts require additional analysis.
<Accumulated contact-mode phase drift> after a duration $T$ is $\Delta\phi=T(\widetilde\omega-\omega)$ at time-grid points. For an original example, let $\omega=1000\,\mathrm{rad/s}$ and $T=0.10\,\mathrm s$. With $h=1\,\mathrm{ms}$, the central-difference <frequency> is about $1047.20\,\mathrm{rad/s}$, producing $4.720\,\mathrm{rad}$ of phase lead. The step satisfies the strict stability inequality, yet the phase lead exceeds half a cycle.
Reducing $h$ to $50\,\mathrm{\mu s}$ gives about $0.01042\,\mathrm{rad}$ of lead over the same duration. At $25\,\mathrm{\mu s}$ it is about $0.002604\,\mathrm{rad}$, approximately one quarter as large, matching the leading second-order error. For a small desired phase budget $\phi_{\max}$, the central-difference expansion suggests
$$
h\lesssim\sqrt{\frac{24\phi_{\max}}{T\omega^3}}.
$$
With $\phi_{\max}=0.01\,\mathrm{rad}$, this estimate is about $49.0\,\mathrm{\mu s}$, much smaller than the $2\,\mathrm{ms}$ stability boundary. Check the exact <frequency> formula when the small-step expansion is not sufficiently accurate.
\Image[/contact-mode-phase-accumulation.png]
{title=Stable contact-mode timesteps can accumulate substantial numerical phase error}
{description=Exact <frequency> errors for central differences and the trapezoidal rule, and the resulting absolute phase errors for a fixed 1000 rad/s mode integrated over 0.1 s. The horizontal line is an illustrative 0.01 rad phase budget. These are analytical linear-mode predictions, not contact measurements.}
In a contact experiment or simulation, a shifted oscillation phase changes when <force> maxima and gap crossings are predicted. Damping, changing tangent <stiffness> and a contact opening can invalidate the fixed-<frequency> model, so the phase estimate is not by itself an event-time error certificate. A proposed convergence check compares <force> histories, dominant <frequency> and event times at successively halved steps, while keeping contact geometry and constitutive parameters fixed. Agreement in maximum penetration alone can miss a phase error that accumulates over many cycles.
Back to article page
