How should I compute the Jacobian for my equations of motion?

Hi there!

I have a set of second-order equations of motion written symbolically in Matlab, put into matrix form using EquationsToMatrix( ). Then I convert these symbolic equations to numerical code files, using matlabFunction( ). With the numerical files, I then wrote six first-order equations of motion, in order to pass them into ode45 for simulations, in a separate simulation file. Now, I have found fixed points in my system, and I want to evaluate their stability. I am supposed to compute the Jacobian partial-derivatives matrix for the right-hand-side functions within the full set of equations of motion. How can I best compute the 6x6 Jacobian, knowing that later I have to find eigenvalues of this Jacobian matrix? Should I work numerically, symbolically, or both? If I work numerically, should I be using the gradient function, or something else?

Thanks in advance,

Answer 1

[Torsten]

According to your description, at some stage of your procedure, the right-hand side of your ODE system is available as a symbolic vector depending on the state variables y. Here, you can easily compute the Jacobian in symbolic form and later on, substitute your fixed points and compute the eigenvalues as shown below.

syms y(t)
eqn = diff(y,2) + y^2 == 4 - diff(y);
[V,S] = odeToVectorField(eqn)

V = 

S = 

% Convert the Y[i] in odeToVectorField output V to yi
V_c = feval(symengine,'evalAt',V,'Y=[y1,y2]')

V_c = 

% Evaluate jacobian with respect to yi
J = jacobian(V_c, sym('y', [1 2]))

J = 

s = symvar(J)

s = 

Jnum = double(subs(J,s(1),2))

Jnum = 2×2

0 1 -4 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

eig(Jnum)

ans = 
  -0.5000 + 1.9365i
  -0.5000 - 1.9365i

[Noob]

Hi Torsten!

My symbolic EOM has 3 second-order equations. In my numerical file, a Matlab function file, I write the corresponding six first-order equations — primarily to pass it to ode45 later, in a separate simulation file. So, I think I need to work numerically, because I need a square 6x6 Jacobian in order to find eigenvalues — my dynamical state has six variables. So, in my symbolic file, I would only have three equations, which wouldn’t be enough. What do you think?

Thanks!

[Torsten]

So, in my symbolic file, I would only have three equations, which wouldn’t be enough. What do you think?

I don't know your code, but using "odeToVectorField" like I did above will give you the six first-order equations from your three second-order equations. That's also the form that is required for ode45 - so I think at some stage you do this transformation.

If you find it easier to do it numerically, you can use

https://uk.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

Function to use should be "jacobianest".

[Noob]

HI Torsten!

Thanks so much for your answer and replies. Let me try to puzzle this out in the coming days.

I had a quick question for you, if you don't mind: Can I switch from symbolic to numerical, back and forth, using sym( ) and double( ), without any consequences, and that Matlab's computations should be trustworthy and precise? Is there something to be aware of, when switching from symoblic to numerical, and vice versa? I'm aware that some Matlab functions are symbolic-only, such as the Jacobian command, so I'm aware of this limitation, at the very least. Anything else to be mindful of?

Thanks again for your help!

[Torsten]

The main source of errors you can encounter is due to the integration of the differential equations with the numerical solver. So don't worry about the switch from symbolic to numeric and vice versa - it's neglectable.