Time stepping in Rigid Body Dynamics

Question

When simulating the interaction of 3D rigid bodies, how do you decide which object to advance its state first? If you advance all bodies at the same time there will be intersection/constraint violation. Maybe the system could then be globally solved for but it seems this is likely to produce weird results.

Answer 1

Indeed, naïvely solving each constraint a single time and sequentially is extremely inaccurate, as doing so generically violates previous constraints and so the first constraint will typically be poorly enforced.

To counter this, there are essentially three ways to update the bodies at each time step:

• Global Constraint Solving: You stack all of the Jacobian matrices for each single constraint into one large $6n\times c$ matrix, stack the generalised velocity vector for each object: $\mathbf v=[\vec v_1\ \vec\omega_1\ \vec v_2\ \vec\omega_2\ ...]^\top$, form the generalised system mass $\mathbf M^{-1}=\mathrm{diag}\left(m_1^{-1}\ m_1^{-1}\ m_1^{-1}\ \mathbf I_1^{-1}\ ...\right)$, and compute the force on every object simultaneously as $\mathbf F=\mathbf J^\top\boldsymbol{\lambda}$ where, employing Baumgarte stabilisation, $\boldsymbol\lambda$ is solved as: $$ (\mathbf{JM}^{-1}\mathbf{J^\top})^{-1}\boldsymbol\lambda=\frac{1}{\Delta t}\boldsymbol \beta-\mathbf J\left(\frac{1}{\Delta t}\mathbf v+\mathbf M^{-1}\mathbf f_\text{ext}\right) \\\sim\mathbf A\boldsymbol{\lambda}=\mathbf b $$ for which you can perform a Cholesky decomposition, Jacobi, or whatever you like. Then you apply the force directly, or opt for the impulse update $\Delta\mathbf v=\mathbf M^{-1}\mathbf J^\top\boldsymbol\lambda$. This will essentially yield perfect results, provided a global solution does indeed exist for all the constraints. However, this method is downright terrible for game physics engines! This is because the size of the stacked Jacobian scales as $\mathcal O(N^2)$, where $N$ is the number of bodies, and so you will have to solve extremely large matrix equations to determine the corrective force, with a time complexity of $\mathcal O(N^{\sim2.8})$. Also, from experience, it can lead to single-frame collapse, where there is no global solution, so it just ends up destroying all the constraints without employing time-consuming mitigation tactics. That said, this method is perfect for pre-built rendered physics, where speed is not a major concern.

• Iterative Constraint Solving: This is somewhat like a mean-field approach. You loop over each constraint, solve it on its own, then feed the results into the sequential constraints. When you iterate multiple times over all of the constraints, the accumulated impulses/forces will converge to the global solution. In other words, for each constraint $C_i$ you compute $\boldsymbol\lambda_i$ multiple times (since it will generically change after the velocity updates from the other constraints). The resulting impulse is accumulated (subject to requisite clamping). Contrary to global solvers, iterative solvers are linear in both space and time: $\mathcal O(N)$. The accuracy will be a bit sketchy in a naïve implementation due to slow convergence to the global solution, but you can improve this through various techniques like warm starting (inter-frame coherence).

• Block Solving: a mix between the previous two, bringing in both of their benefits. Note that in a lot of cases, you actually do want to solve multiple constraints simultaneously for stability reasons - the best example I can think of is a resting inequality constraint, where solving them sequentially will lead to unwanted oscillation of the top rigidbody. You identify "blocks" of objects that you want to solve together, and do so for each block is multiple iterative runs. This is very useful for accommodating blocks of sleeping rigid bodies too.

Answer 2

You need to step them all at the same time.

You might

• Identify all such collision/constraint violations
• Do those steps again differently. Figure out when the collision occurred and step just to that point. Figure out the rebound. Finish the step in a new direction.

It would help to make time steps small enough that you don't have double bounces in a single step.