Equality Constraints

Impulse and Force Based Constraints

\mathbf{a} = \begin{bmatrix}a_x\\ a_y\\ \alpha\end{bmatrix} = \mathbf{M}^{-1}\mathbf{F} = \begin{bmatrix}\frac{1}{m} & 0 & 0\\0 & \frac{1}{m} & 0\\0 & 0 & \frac{1}{I}\end{bmatrix}\begin{bmatrix}F_x\\ F_y\\ \tau\end{bmatrix} = \begin{bmatrix}\frac{1}{m} & 0 & 0\\0 & \frac{1}{m} & 0\\0 & 0 & \frac{1}{I}\end{bmatrix}\begin{bmatrix}F_x\\ F_y\\ (\mathbf{p}-\mathbf{c})\times \begin{bmatrix}F_x\\ F_y\end{bmatrix}\end{bmatrix}
Velocities are changed by impulses via \mathbf{v}:=\mathbf{v}+\Delta\mathbf{v}, with
\Delta\mathbf{v}=\mathbf{M}^{-1}\mathbf{F}\Delta t=\mathbf{M}^{-1}\mathbf{P}=\begin{bmatrix}\frac{1}{m} & 0 & 0\\0 & \frac{1}{m} & 0\\0 & 0 & \frac{1}{I}\end{bmatrix}\begin{bmatrix}P_x\\ P_y\\ P_\alpha\end{bmatrix}=\begin{bmatrix}\frac{1}{m} & 0 & 0\\0 & \frac{1}{m} & 0\\0 & 0 & \frac{1}{I}\end{bmatrix}\begin{bmatrix}P_x\\ P_y\\ (\mathbf{p}-\mathbf{c})\times \begin{bmatrix}P_x\\ P_y\end{bmatrix}\end{bmatrix}

Computing Constraint Impulses

1. To ensure the constraint,
\dot{C}=\mathbf{J}\mathbf{v}+b=0
2. Constraint forces or impulses may not add energy to the system, i.e. do work - force times distance: \mathbf{F}^T\Delta\mathbf{x}=\mathbf{F}^T\Delta t\mathbf{v}=\mathbf{P}^T\mathbf{v}=0. Since \mathbf{J}\mathbf{v}=0, the impulse we search must be a multiple of the transposed Jacobian:
\mathbf{P}=\lambda\mathbf{J}^T
NB, some tutorials use the bead-on-a-wire example and state that the constraint force or impulse stands orthogonally on the velocity. This may lead to a misconception; I wouldn't attach much geometrical meaning to \mathbf{P}^T\mathbf{v}=0, since both are accumulated vectors. For a joint constraint, the individual impulses roughly point from one point to the other.
3. As discussed further above, the impulse is applied to the previously predicted, tentative velocity \overline\mathbf{v} via
\mathbf{v} = \overline{\mathbf{v}} + \mathbf{M}^{-1}\mathbf{P}

Example: Complete Joint Constraint Derivation

Constraint Solver

for i = 1 to nIterations
    for c in constraints
        P = computeImpulses(c);
        c.bodies.applyImpulses(P);

Example

tStop:

t=0s

Notes

Note, I did not add parameters to this simulation; frankly, because using links with less mass leads to jitter and I didn't bother taking care of it. I simply went to 1/4th of the previous time step size - 4ms. I suppose that rope simulations with a larger number of links always lead to problems; if one uses not sequential impulses but springs to simulate it, one gets a stiff system of equations. There's probably a way to take better care of it though; but keep in mind a rope is a rather extreme application of joints because of the many couplings. Other, similarly complex applications of joint chains like ragdolls are treated specifically as far as I know.