You can also do this in Box2D (press P to toggle position correction). Of course we can't guarantee the model will be assembled in the correct kinematic branch. You can see this in demo 9 where sometimes the links will get twisted into a knot.
I just tried to start Box2D but it doesn't work with my glut32.dll. Doesn't matter.
I try to explain the difference between your method and mine. Please correct me, if I say something of wrong.
Let us assume you want to simulate two bodies connected by a ball joint. Then the position constraint is a-b=0 where a and b are the joint points.
Erin's method:
To simulate this constraint you differentiate the constraint equation with respect to time in order to get a constraint for the velocities. So in fact the position constraint is satisfied by keeping the velocities of the joint points equal. Then you compute the Jacobian J and solve the equation
M dV/dt = J^T lambda + Fext
for lambda. The problem is that the two joint points will drift apart due to numerical errors. Since the position of these points is not regarded in the velocity constraint it can't be corrected directly. So a stabilization method like the one from Baumgarte is required.
My method:
I don't differentiate the position constraint. The first step of my method is that a prediction for the positions of both joint points is made by integrating the equation
r' = omega x r
where r is the vector from the center of mass c to a joint point. The new positions at time t+h are
a(t+h) = r_a(t+h) + c(t+h)
b(t+h) = r_b(t+h) + c(t+h)
This prediction determines the positions of the joint points, if the bodies have a ballistic motion (without regarding the constraint). The distance at time t+h should be zero. So the constraint I try to satisfy looks like this:
d(t+h) = b(t+h)-a(t+h) = 0
Since the motion of the joint points is almost linear in a small time step, I use a linearization. So I determine an impulse that changes the velocities by d(t+h)/h by solving equation
K(t) * p = d(t+h)/h
where (here J is the inertia tensor not the Jacobian)
K = (1/m1+1/m2) I3 - (r_a* J1^-1 r_a* + r_b* J2^-1 r_b*)
Since K is symmetric, positive definite, constant for time t and regular (proof by Mirtich), the equation can be solved by inverting K. The resulting impulse p is applied to the bodies in opposite directions at time t (not t+h). Due to the velocity change at time t the distance d(t+h) gets smaller. In the next iteration the predicted value d(t+h) has to be computed again. Since K is constant for time t, the matrix has not to be inverted again. So in the second iteration the prediction and a simple matrix multiplication is all we need to compute the impulse.
One advantage of my method is that no additional stabilization method like Baumgarte is required, since the future distance d(t+h) is always zero after the time step.
Velocity constraint can be handled in an easier way since no prediction is required.
I hope that helps to understand the difference between the methods.
Jan