I was reading this article:

https://www.toptal.com/game/video-game-physics-part-iii-constrained-rigid-body-simulation

and it is missing the derivation of a Jacobian.

Can somebody please do the derivation of the Jacobian with two equality constraints for double pendulum?

I want to simulate double pendulum by using force-based approach with system solver, not the iterative approach.

The author gives an explanation of the Jacobian. The Jacobian rows are the transposed gradients of the scalar constraint functions. In practice the Jacobian is usually found 'by inspection' using that the time derivative of the constraint function has a special form, dC/dt = J*v 

E.g. say you have a simple distance constraint between two particles 

C =| p2 - p1|

dC/dt =  (p2 - p1)/ |p2 - p1| * (v2 - v1) = n * (v2 - v1) , where n = (p2 - p1) / | p2 - p1|

J = ( -n^t n^t )

Going through the time derivative is often way easier than building partial derivatives/gradients...

I don't know how to do it.

Can you please derive two equality constraints for double pendulum?

You mean like this? http://scienceworld.wolfram.com/physics/DoublePendulum.html

This is reduced coordinates using Euler Lagrange equation.

I did. It is the same as I wrote above.

C1 = |p2 - p1| - L

C2 = |p3 - p2| - L

J1 = ( -n1^t | n1^t | 0^t ) where n1 = (p2 - p1) / |p2 - p1|

J2 = ( 0^t | -n2^t | n2^t ) where n2 = (p3 - p2) / |p3 - p2|

On 6/10/2019 at 6:20 PM, Randy Gaul said:

You mean like this? http://scienceworld.wolfram.com/physics/DoublePendulum.html

This is reduced coordinates using Euler Lagrange equation.

I think it should be done without using angles, just positions.

Just like my example article solves single pendulum just with positions.

BTW what kind of math is this (in my example article)?

Is it vector calculus or analytical mechanics or what, can you please tell me?

It's a mix of vector calc, classical mechanics, and linear algebra. It all comes from Erin's online resources here: http://box2d.org/downloads

I am incline to doubt that only video games use these kind of a simulations.

I meant the article you're reading is largely copied from Erin's stuff.

I figured it out partially, but it doesn't work for masses other than 1.0

Can you please tell me how to fix this?

        Mx J = new Mx(2, 4);
        Mx Jdot = new Mx(2, 4);
        Mx M = new Mx(4);
        vector Fext = new vector(4);
        vector qdot = new vector(4);
        vector lambda;


                M.Zero();
                M.data[0, 0] = 1.0 / tp.mass;
                M.data[1, 1] = 1.0 / tp.mass;
                M.data[2, 2] = 1.0 / tp3.mass;
                M.data[3, 3] = 1.0 / tp3.mass;

                Fext.zero();
                Fext.data[0] = 0;
                Fext.data[1] = body.g;
                Fext.data[2] = 0;
                Fext.data[3] = body.g;

                qdot.zero();
                qdot.data[0] = tp.velocity.x;
                qdot.data[1] = tp.velocity.y;
                qdot.data[2] = tp3.velocity.x;
                qdot.data[3] = tp3.velocity.y;

                Jdot.Zero();
                Jdot.data[0, 0] = tp.velocity.x;
                Jdot.data[0, 1] = tp.velocity.y;
                Jdot.data[1, 0] = tp.velocity.x - tp3.velocity.x;
                Jdot.data[1, 1] = tp.velocity.y - tp3.velocity.y;
                Jdot.data[1, 2] = tp3.velocity.x - tp.velocity.x;
                Jdot.data[1, 3] = tp3.velocity.y - tp.velocity.y;

                J.Zero();
                J.data[0, 0] = tp.position.x;
                J.data[0, 1] = tp.position.y;
                J.data[1, 0] = tp.position.x - tp3.position.x;
                J.data[1, 1] = tp.position.y - tp3.position.y;
                J.data[1, 2] = tp3.position.x - tp.position.x;
                J.data[1, 3] = tp3.position.y - tp.position.y;

                lambda = Mx.Invert(J * M * !J) * (-Jdot * qdot - J * M * Fext);


                var fvec = !J * lambda;