I am trying to find the Lyapunov exponents at a fixed point as given in the paper here(page 14). The code that I have used till now is given below:

K1 = 6.90; K2 = 16.05; K3 = 9.65; K4 = 3.27; K5 = 6.55;
\[Omega]sq[0] = -0.923; \[Omega]sq[1] = 6.478;
lagrangian = 
  Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] + 
   K1 c[0][t]^3 + K2 c[0][t] c[1][t]^2 + K3 c[0][t] c[0]'[t]^2 + 
   K4 c[0][t] c[1]'[t]^2 + K5 c[0]'[t] c[1][t] c[1]'[t];
c[0][t_] := OverTilde[c][0][t] + \[Alpha]1*OverTilde[c][0][t]^2 + \[Alpha]2*OverTilde[c][1][t]^2; 
c[1][t_] := OverTilde[c][1][t] + \[Alpha]3*OverTilde[c][0][t]*OverTilde[c][1][t]; 
\[Alpha]1 = -2; \[Alpha]2 = -0.5; \[Alpha]3 = -1;
n = Expand[lagrangian];
vars = {OverTilde[c][0][t], OverTilde[c][1][t], 
lagrangian = 
  Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1;
momentum[n_] := D[lagrangian, Derivative[1][OverTilde[c][n]][t]]
hamiltonian = Expand[Sum[momentum[n]*Derivative[1][OverTilde[c][n]][t], {n, {0, 1}}] - lagrangian]; 
eulerLagrange[lagrangian_, vars_, dvars_] := 
       Thread[Table[D[D[lagrangian, dvar], t], {dvar, dvars}] - Table[D[lagrangian, var], {var, vars}] == 
         ConstantArray[0, Length[vars]]]; 
    equationsOfMotion = eulerLagrange[lagrangian, {OverTilde[c][0][t], OverTilde[c][1][t]}, 
        {Derivative[1][OverTilde[c][0]][t], Derivative[1][OverTilde[c][1]][t]}]

1 Answer 1


After having a look at the code you have and the paper you linked, it seems to me that to get to your goal, you need the following steps (This is too long for a comment):

  1. Transform your differential equations to a first order system of differential equations, so your equations are in the form dx/dt = F(x).

  2. Once you have the vector F(x), you can calculate its Jacobian in mathematica with the function: jacobianMatrix[f_List, x_List] := Outer[D, f, x], e.g.
    j=jacobianMatrix[f, {c1,c2,c3,c4}] (or however you want to name your variables)

  3. For the dynamics of a local fixed point, the Lyapunov exponents are just the eigenvalues of the fixed point. You will first need to find the fixed point location via something like fp = NSolve[f==0,{c1,c2,c3,c4}]. Then plug the fixed point location into your jacobian matrix j/.fp and finally calculate its eigenvalues with eigenvalues = Eigenvalues[j/.fp].


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.