Joint's Jacobian Matrix is derived from Energy Function[1], for joint solver I followed Erin Catto’s iterative method called Projected Gauss Seidel (PGS) [2] with some modification, call it Projected SOR with prediction and correction to make it more stable for stiff joint.
Erin Catto’s Sequential Impulse[3], which is more suitable for contact, is mathematically equivalent to PGS, and is not implemented current.
Reference
[1] David Barraf. Physic Based Modeling Course Note
[2] Erin Catto. Iterative Dynamic
[3] Erin Catto. GDC2007 Modeling and Solving Constraints