Chain simulation (sequential vs simultaneous solver)

Intro

This is a quick exercise to learn how constraints work in a physics simulator. The best way to learn something is to teach it, so here I go.

This physics simulator is heavily based on Erin Catto's GDC2009 talk.

Simulation loop

The simulation loop looks like this:

• Integration: Computes tentative velocities: $$v_{tentative} = v_i + a \cdot t $$
• Constraints enforcement: Corrects tentative velocities so they meet the constraints, for this we need to figure out how the constraint looks like in order to compute the Jacobian. The Jacobian will help us compute the velocity correction. The Jacobian is a matrix that we'll need to invert. Depending on the solving method this J will be bigger or smaller.
  
  These are the 2 methods we are using here:
• Simultaneously: all the 'distance constraints' actuate at once, in this case the Jacobian is a big matrix, see the 'Simultaneous solver Jacobian computation', this is slow but precise
• Sequentially: applying 1 'distance constraint' each link of the chain and repeating this a few times. In this case the Jacobian is a tiny matrix, this is fast but good for games :)

In every case we will have to:
• Compute the Jacobian
• Uses a velocity bias to prevent drift/stretching (Baumgarte stabilization): $$b = \frac{\beta}{t} \cdot C$$
• Compute the inverse of $JM^{-1}J^T$ (this is why we want smaller jacobians).
• Compute lagrange multipliers (anything that multiplies a derivative that points in the direction we want is called a lagrange multiplier, this is just to scare people).
$$\lambda = -\frac{J \cdot v_{tentative} + b}{JM^{-1}J^T}$$
• And the velocity gets corrected(so it meets the constraints) by applying:
$$v_f = v_{tentative} + \lambda J^T$$
• Rendering: Draws particles and links.

Simultaneous solver Jacobian computation

For a 4 element chain the constraints are: $$\begin{align}C_1(p1,p2,p3,p4)=& .5 \cdot |p_1-hook|^2 -d^2 =0 \\ C_2(p1,p2,p3,p4)=& .5 \cdot |p_1-p_2|^2 -d^2 =0 \\ C_3(p1,p2,p3,p4)=& .5 \cdot |p_2-p_3|^2 -d^2 =0 \\ C_4(p1,p2,p3,p4)=& .5 \cdot |p_3-p_4|^2 -d^2 =0 \\ \end{align}$$ Where:

• $hook$ is the point were the chain hangs from.
• $p_i$ are the position of each particle.
• $d$ is the distance of each link.

For each $C_i$ we need to compute $\frac{dC_i}{dt}$, using the total derivative and the chain rule we have that: $$\frac{dC_i}{dt}=\frac{\partial C_i}{\partial p_1}\cdot\frac{dp_1}{dt} + \frac{\partial C_i}{\partial p_3}\cdot\frac{dp_3}{dt} + \frac{\partial C_i}{\partial p_4}\cdot\frac{dp_4}{dt}$$ Where: $$\frac{\partial C_i}{\partial p_1} = \frac{\partial (p_i-p_{i-1})^2 }{\partial p_i} = 2 \cdot (p_i-p_{i-1}) \qquad\text{and}\qquad \frac{dp_i}{dt} = v_i$$ This produces the following: $$\begin{array}{rr} \frac{dC_1}{dt}=&(p_1-hook)\cdot{v_1}& + 0\cdot v_2& + 0\cdot v_3& + 0\cdot v_4& \\ \frac{dC_2}{dt}=& (p_1-p_2)\cdot{v_1}& - (p_1-p_2)\cdot v_2& + 0\cdot v_3& + 0\cdot v_4& \\ \frac{dC_3}{dt}=& 0\cdot v_1& + (p_2-p_3)\cdot v_2& - (p_2-p_3)\cdot v_3& + 0\cdot v_4& \\ \frac{dC_4}{dt}=& 0\cdot v_1& + 0\cdot v_2& + (p_3-p_4)\cdot v_3& -(p_3-p_4)\cdot v_4& \\ \end{array}$$ Which in a matrix form it looks like this: $$\begin{pmatrix}p_1 & 0 & 0 & 0 \\ (p_1-p_2) & -(p_1-p_2) & 0 & 0 \\ 0 & (p_2-p_3) & -(p_2-p_3) & 0 \\ 0 & 0 & (p_3-p_4) & -(p_3-p_4) \\ \end{pmatrix} \cdot \begin{pmatrix} v_1 \\ v_2 \\ v_3 \\ v_4\end{pmatrix}$$ The matrix on the left is the Jacobian.

Contact/Questions:

<my_github_account_username>$@gmail.com$.