= Contact-mode timestep stability
<Contact-mode timestep stability> concerns whether a numerical integrator introduces growing errors when advancing the fast dynamics of an active penalty contact. It depends on the integration rule and the local modal <frequencies>, not just the size of the modeled overlap. The primary compliant-contact treatment identifies stiff contact springs as a source of numerically stiff differential equations.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, Compliant Contact Models, stiff-spring and numerical integration discussion]{ref}
Freeze the active contacts and linearize about a compressed equilibrium. An undamped mode has $\ddot x+\omega^2x=0$. For step $h$, put $z=h\omega$. The <Forward Euler method> updates displacement and velocity with matrix
$$
G_E=\begin{pmatrix}1&h\\-h\omega^2&1\end{pmatrix},
$$
whose eigenvalues are $1\pm iz$. Their magnitudes are $\sqrt{1+z^2}>1$ whenever $z>0$. Thus the method has no nonzero stable step for indefinite integration of this undamped oscillator, although smaller steps reduce its error on a fixed finite interval.
The central-difference recurrence instead gives
$$
x_{n+1}-2x_n+x_{n-1}=-z^2x_n,
\qquad r^2-2(1-z^2/2)r+1=0.
$$
For $0<z<2$ its two distinct amplification roots lie on the unit circle. At $z=2$ the repeated root is $-1$, and solutions can include $n(-1)^n$, so the strict bound matters. A collection of frozen undamped modes therefore requires $h\omega_{\max}<2$ for this oscillatory stability condition. In <coupled normal-contact penalty modes>, $\omega_{\max}$ comes from the coupled system rather than an arbitrary total-<mass> estimate.
Stability still allows substantial <contact-mode numerical phase error>. For central differences the phase per step is $\theta=2\arcsin(z/2)$, giving numerical angular <frequency> $\theta/h$. At $\omega=1000\,\mathrm{rad/s}$ and $h=1\,\mathrm{ms}$, $z=1$ is stable but the numerical <frequency> is about $1047.2\,\mathrm{rad/s}$, $4.72\%$ high. A timing-sensitive contact event can consequently drift even when oscillation amplitudes do not grow.
Implicit schemes illustrate another tradeoff. For the <Backward Euler method>, the amplification factors are $(1\mp iz)^{-1}$ and have magnitude $1/\sqrt{1+z^2}$. At $z=1$, each step reduces the oscillatory amplitude by about $29.3\%$, despite the continuous undamped system conserving energy. The trapezoidal rule instead has factors
$$
r_\pm=\frac{1\pm iz/2}{1\mp iz/2},\qquad |r_\pm|=1,
$$
with phase $2\arctan(z/2)$. The same example gives $927.3\,\mathrm{rad/s}$, $7.27\%$ low. Unconditional stability for this linear model does not imply correct timing at a large step.
A practical error check therefore compares both amplitude and phase while reducing $h$. For instance, a small target dimensionless step $h\omega_{\max}=0.1$ gives much smaller phase error than one, but the acceptable value depends on the observable and time span; it is not a universal accuracy certificate.
The frozen-mode argument has a defined scope. Contact opening changes the active <stiffness>, nonlinear <force> laws change tangent <stiffness>, and damping introduces additional eigenvalue scales. An integrator can satisfy the local oscillator bound while missing an intervening collision or contact release. A proposed test starts with a small compressed-contact perturbation, checks <frequencies> and energy behavior, then repeats across a known contact-opening event with decreasing steps. Record <force> histories and event times as well as penetration, so numerical damping or phase drift is not mistaken for a material response.
Back to article page
