the mass matrix and Jacobian matrix in ode15s

Hi there, I have a set of odes in the form

where M and K are n*n matrixes q is a n*1 vector containing the generalised coordinates, and F is the n*1 generalised force vector.

I am using ode15s to solve this problem. Here I want to ask, are M and K the mass matrix and Jacobian matrix, respectively that should be input in odeset?

Any help is geatly appreciated.

Cheers

Answer 1

[Sam Chak]

Hi @Tony Cheng

is the state-dependent mass matrix. Therefore, you can specify the mass matrix using the 'Mass' option of odeset(). However, is not a Jacobian matrix. By the way, please check if the highest degree of the ODE is or .

[Tony Cheng]

Hi Sam, thx for your useful reply! The original odes are in the form of

and the highest degree of the ODE is actually .

To make them suitable for the application of ode solver in Matlab, I rewrote them to the form of

Then, in this problem, is K1 the Jacobian matrix? If not, what the Jacobian matrix is？ Many thanks!

[Torsten]

The Jacobian J is

J = d/dq (M^(-1)*(F-K*q))

Thus if K, M and F do not depend on q,

J = -M^(-1)*K

[Sam Chak]

Hi @Tony Cheng

Here is an example where the state differential equations are given.

The Jacobian matrix, 𝒥, is determined by getting all first-order partial derivatives of the multivariate functions:

.

Check this simple MATLAB script.

syms f g x y
f =   y;
g = - 0.006*y - 10*sin(x);
J = jacobian([f, g], [x, y])

J = 

[Tony Cheng]

Hi Sam,

From your example and Torsten's answer, I find that the Jacobian matrix in odes in the form of

is determined by

. So is it rite? Many thx!

[Sam Chak]

Hi @Tony Cheng,

Try obtaining the state-space model.

Equations of motion

Define the state vector

The state-space becomes

In compact form (assuming that )

Now you can compute the Jacobian matrix.

[Torsten]

@Sam Chak

Maybe you could replace q and qdot on the right-hand side by x1 and x2.

If you set then x = [x1;x2], the Jacobian with respect to x of the right-hand side should be more obvious for the OP.

@Tony Cheng

In short:

The Jacobian of your system is [0, I ; -M_1^(-1)*K_1, -M_1^(-1)*C_1] if your original equation reads M_1*q'' = -K_1*q - C_1*q' + F1 and you substitute [x(1),...,x(n)] = q, [x(n+1),...,x(2*n)] = q'.

[Tony Cheng]

Hi Sam, thx so much for your valuable help! I get the definition of Jacobian matrix now!

[Bruno Luong]

Be careful the formula you get [0, I ; -M_1^(-1)*K_1, -M_1^(-1)*C_1]

only valid for M and K that do not depend on q and qdot.

Otherwise you have to derive this term wrt q and qdot and multiply to [q; qdot] then add to the above expression. It will be messy to write it down with all the terms.

[Tony Cheng]

Hi Bruno, many thx for your reminding.

To make my problem suitable for the ode15s solver, the original problem

(1)

is reformated as

(2)

let , then Eq.(1) can be rewritten as

then the Jacobian matrix is Jacobian(A-BQ, Q). Am I right?

[Torsten]

You have to substitute q = x(1:n) and qdot = x(n+1:2*n) and rewrite your system as

x(1:n)_dot = x(n+1:2*n)
M(x(1:n))*x(n+1:2*n)_dot = -C(x(1:n),x(n+1:2*n))*x(n+1:2*n) - K(x(1:n))*x(1:n) + F(x(1:n))

or

x(1:n)_dot = x(n+1:2*n)
x(n+1:2*n)_dot = (M(x(1:n)))^(-1)*(-C(x(1:n),x(n+1:2*n))*x(n+1:2*n) - K(x(1:n))*x(1:n) + F(x(1:n)))

This system is of the form

x' = F(t,x)

with F being a vector of 2*n functions f(1),...,f(n),f(n+1),...,f(2*n) in 2*n unknowns x(1),...x(n),x(n+1),...,x(2*n).

The (i,j)-th element of the Jacobian J of F is as usual J_ij = df_i/dx_j.

From your notation it seems that M,C,K and F are not constant, but depend on q (and C even on qdot). This makes the computation of the Jacobian very messy, as @Bruno Luong noted correctly.

But fortunately, the stiff ODE integrators will compute the Jacobian numerically and you don't have to care about it.

If in your notation from above M, C and K are really constant and don't depend on q and qdot as you wrote before, the Jacobian is -B.

[Bruno Luong]

@Tony Cheng "then the Jacobian matrix is Jacobian(A-BQ, Q). Am I right?"

No it's -d(B(Q)*Q)/dQ

[Torsten]

"A" might also depend on Q ...

[Bruno Luong]

Oh yeah correct

d(A(Q)-B(Q)*Q) / dQ

Answer 2

[Sam Chak]

Hi @Tony Cheng

Following @Torsten's advice, perhaps showing you an example of obtaining the Jacobian matrix for a 2-DOF rigid robotic manipulator would be the best way for you to learn.

Rearranging the equation to obtain:

syms x1 x2 x3 x4 tau1 tau2

r1 = 1.0;   % length of link 1
r2 = 0.8;   % length of link 2
m1 = 1.0;   %   mass of link 1
m2 = 1.5;   %   mass of link 2

% State vector
x   = [x1;  % q1
       x2;  % q2
       x3;  % q1dot
       x4]; % q2dot

% Mass matrix
M11 = (m1 + m2)*r1^2 + m2*r2^2 + 2*m2*r1*r2*cos(x2);
M12 = m2*r2^2 + m2*r1*r2*cos(x2);
M21 = M12;
M22 = m2*r2^2;
M   = [M11 M12;
       M21 M22]

M = 

% Coriolis matrix
C12 = m2*r1*sin(x2);
C   = [-C12*x4 -C12*(x3+x4);
        C12*x1  0]

C = 

% Gravity
G1  = (m1 + m2)*r1*cos(x2) + m2*r2*cos(x1 + x2);
G2  = m2*r2*cos(x1 + x2);
G   = [G1;
       G2]

G = 

% Control torque vector
T   = [tau1; 
       tau2];

% Dynamics in state-space form, xdot = F(x, τ)
ddq = inv(M)*(T - C*[x3; x4] - G);
f1  = x3;
f2  = x4;
f3  = ddq(1);
f4  = ddq(2);
F   = [f1;
       f2;
       f3;
       f4];

% Jacobian matrix
J = jacobian(F, x)

J = 

[Bruno Luong]

I made some calculation of the ODE Jacobian and I put in the attached pdf.

Warning : there might be some mistake about it because I do not make any check and it's all about pencil+paper.

EDIT: I'm pretty much certain I get right the formula now.

[Tony Cheng]

Hi Sam Chak, Bruno Luong, and Torsten, thanks so so much for your great patience, help and time in giving me so many useful answers for this problem! I think the ode15s solver can be much much better used in my project with a right Jacobian in the odeset!

[Sam Chak]

Thanks @Bruno Luong and @Torsten for showing the analytical calculations of the Jacobian matrix. A good learning experience for @Tony Cheng and me.

[Walter Roberson]

@Tony Cheng commented to @Bruno Luong about the same time as the "hi" posting:

Hi Bruno, the pdf file is very useful, caz the M(q), C(\dot{q},q), K(q), and F(q) in my porject is very dense, thus, the final details can contribute to a simpler symbolic calculation in Matlab.

[Bruno Luong]

I correct the pdf file of Jacobian formula (one of the kronecker is reversed)

[Sam Chak]

@Bruno Luong, Thanks for the update on the Jacobian matrix for the Euler–Lagrange equations of motion in matrix form.

[Walter Roberson]

On August 14, @Tony Cheng commented to @Torsten

Hi Torsten, thx so much for your valuable help! I get the definition of Jacobian matrix now!

Answer 3

[Bruno Luong]

This script tests the formula with dummy test configuration by comparing with the Jacobian obtained from finite difference:

runtest

Jdiff = 6×6

0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 1.9007 -1.0509 -0.9406 -0.0083 -0.0099 -0.0155 -0.4763 3.2214 -0.6506 -0.0057 -0.0068 -0.0107 -0.5330 -0.8133 2.7854 -0.0064 -0.0076 -0.0120

J = 6×6

0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 1.9007 -1.0509 -0.9406 -0.0083 -0.0099 -0.0155 -0.4763 3.2214 -0.6506 -0.0057 -0.0068 -0.0107 -0.5330 -0.8133 2.7854 -0.0064 -0.0076 -0.0120

function runtest

n = 3;
q = rand(n,1);
qd = rand(n,1);

hq = 1e-6;
hqd = 1e-6;
[y, M, K, C, F, dMdq, dKdq, dCdq, dCdqd, dFdq] = odefun(q, qd);

% Compute the Jacobian by finite difference
n = length(q);
Jdiff = zeros(2*n);
for j=1:n
    ej = accumarray(j,1,[n 1]);
    yq = odefun(q+hq*ej, qd);
    yqd = odefun(q, qd+hqd*ej);
    Jdiff(:,j) = (yq-y)/hq;
    Jdiff(:,n+j) = (yqd-y)/hqd;
end
Jdiff

% Compute the Jacobian by the formula in the pdf paper
I = eye(n);
z = zeros(n);
qdd = -(M \ (K*q + C*qd - F));
Dq = odot(dMdq,qdd) + odot(dKdq,q) + odot(dCdq,qd) - diag(dFdq);
J = -[z, -I;
      M\(K+Dq), M\(C+odot(dCdqd,qd))];
J
end

function AB = odot(A,B)
% Implement
% odot = @(A,B) A*kron(eye(width(A)/height(B)), B);
A = reshape(A, size(A,1), size(B,1), []);
AB = pagemtimes(A, B);
AB = reshape(AB, size(A,1), []);
end

%% Function that simulate the ODE state y = -dQ/dt = -[qdot; qdotdot]
function [y, M, K, C, F, dMdq, dKdq, dCdq, dCdqd, dFdq]  = odefun(q, qd)
[M, dMdq] = Mfun(q);
[K, dKdq] = Kfun(q);
[C, dCdq, dCdqd] = Cfun(q, qd);
[F, dFdq] = Ffun(q);
I = eye(size(M));
z = zeros(size(M));
Mtilde = [I, z;
          z  M];
Ktilde = [z -I;
          K  C];
Q = [q; qd];
Ftilde = [z(:,1); F];
y = Mtilde \ (Ftilde-Ktilde*Q);
end

function A = setdiag(v)
n = length(v);
A = zeros(n,n,n);
for k=1:n
    A(k,k,k) = v(k);
end
A = reshape(A, n, n*n);
end

%% functions that generate Mass, Stiffness, Coriolis and forcing
function [M, dMdq] = Mfun(q)
M = diag(1 + q.^2);
dMdq = setdiag(2.*q);

a = 0.1;
s = a*sum(q.^2);
M = M + s;
n = length(q);
dsdq = 2*a*q;
dsdq = repelem(reshape(dsdq,1,n),n,n);
dMdq = dMdq + dsdq;

end

%%
function [K, dKdq] = Kfun(q)
K = diag(1 + sqrt(q));
dKdq = setdiag((1/2)./sqrt(q));

a = 0.1;
s = a*sqrt(sum(q.^2));
K = K + s;
n = length(q);
dsdq = a*q./sqrt(sum(q.^2));
dsdq = repelem(reshape(dsdq,1,n),n,n);
dKdq = dKdq + dsdq;

end

%%
function [C, dCdq, dCdqd] = Cfun(q, qd)
C = 0*diag(sin(q).*cos(qd));
dCdq = 0*setdiag(cos(q).*cos(qd));
dCdqd = 0*setdiag(-sin(q).*sin(qd));

a = 0.1;
b = 1/4;
c = b*sum(qd.^2);
s = a*sin(c);
C = C + s;
n = length(q);
dsdq = 2*a*b*cos(c)*qd;
dsdq = repelem(reshape(dsdq,1,n),n,n);
dCdqd = dCdqd + dsdq;

end

%%
function [F, dFdq] = Ffun(q)
F = exp(2*q);
dFdq = 2*exp(2*q);
end

[Tony Cheng]

Hi Bruno, the numerical results of Jdiff and J are identical with four decimals in this case. In my work, I am not sure the results will be the same with or without offering the Jacobian function in odeset. I will try to find whether this term has influences on the final results from ode15s.

BTW, many thanks for ur codes. Really amazing!

[Bruno Luong]

Normally the analytic formula is more accurate. In practice I do not think it makes a big difference between analytic formula and finite difference, unless if your ode system has sate vector with large dimension.

In robotics it has about 6 joints no?

It is just fun to derive the Jacobian for such system. Such formula might help to analyze where and when the system becomes singular. The final formula that has the same pattern as the ode, even the intermediate calculation seems messy and technical not easy.

But I'm curious to know if it helps ode15s to set 'Jacobian' in odeset.

[Tony Cheng]

Hi Bruno, yeah, it is techanically very difficult to derive the analytical expression of the Jacobian matrix. However, in the ode15 solver, it is said that a Jacobian matrix in odeset is very important for stiff odes. I will try to see if it helps. Hope it can play a critical role in the numerical computation. The size of q in my project is 16. I have no idea that the system is large or moderately large in Matlab. Besides, the matrixes M, C and K are very dense. Many of their entries are very lengthy and complicated.

.