I am trying to reproduce the convergence plot of the four Lyapunov exponents for a string from this paper (page 12, figure 7).

The code that I have used till now to find the equations 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]}]
  • $\begingroup$ If you can convert your second-order system to a system of first-order equations, you could use the LyapunovExponents function here. I tried, but couldn't even get the system to NDSolve on it own, due to stiffness. Did you manage to numerically solve your system? What are the initial conditions? $\endgroup$
    – Chris K
    May 2 at 20:06
  • $\begingroup$ @ChrisK, at first I tried solving my system with Marco Sandri's package for Lyapunov Exponents, but the code kept on running indefinitely without any output. When I tried your function, after converting my system to a first-order one, I am getting lots of errors and warning messages. The initial conditions for the problem are specified in figure 7 of this paper: arxiv.org/pdf/2111.09441.pdf (on page 12) $\endgroup$
    – codebpr
    May 3 at 8:06
  • $\begingroup$ Can you get even NDSolve to work on the model? Before calculating the Lyapunov exponents, you should be able to solve the dynamics. $\endgroup$
    – Chris K
    May 3 at 15:29
  • $\begingroup$ Actually, I used the above equations to plot a Poincare section using NDSolve and WhenEvent, hence I assumed it to be working. $\endgroup$
    – codebpr
    May 4 at 6:25

1 Answer 1


I uncovered an unexpected little bug in LyapunovExponents that caused those errors. If you use the updated version from here, it now seems to work.

Wrangle into first-order form:

eqns = Flatten@Join[
      equationsOfMotion, {Derivative[2][OverTilde[c][0]][t], 
       Derivative[2][OverTilde[c][1]][t]}] /. {
      Derivative[2][OverTilde[c][0]][t] -> Derivative[1][dc[0]][t], 
      Derivative[2][OverTilde[c][1]][t] -> Derivative[1][dc[1]][t],
      Derivative[1][OverTilde[c][0]][t] -> dc[0][t], 
      Derivative[1][OverTilde[c][1]][t] -> dc[1][t]} /. (lhs_ -> 
       rhs_) -> (lhs == rhs),
   {OverTilde[c][0]'[t] == dc[0][t], OverTilde[c][1]'[t] == dc[1][t]}];

Warm up system to get good ICs for LyapunovExponents (note: I guessed dc[1][0]==0 since the paper didn't give it explicitly):

tmax = 100;
sol = NDSolve[{eqns, {OverTilde[c][0][0] == -0.0002, 
      OverTilde[c][1][0] == 0.0011, dc[0][0] == 0, 
      dc[1][0] == 0}}, {OverTilde[c][0], OverTilde[c][1], dc[0], 
     dc[1]}, {t, 0, tmax}][[1]];
Plot[Evaluate[{OverTilde[c][0][t], OverTilde[c][1][t], dc[0][t], 
    dc[1][t]} /. sol], {t, 0, tmax}, PlotRange -> All]
ics = {OverTilde[c][0] -> (OverTilde[c][0][tmax] /. sol), 
  OverTilde[c][1] -> (OverTilde[c][1][tmax] /. sol), 
  dc[0] -> (dc[0][tmax] /. sol), dc[1] -> (dc[1][tmax] /. sol)};

enter image description here

Calculate the exponents:

LyapunovExponents[eqns, ics, ShowPlot -> True, PlotExponents -> 4, 
 PlotOpts -> {AxesLabel -> {"step", "\[Lambda]"}, 
   PlotRange -> {-0.006, 0.006}, GridLines -> Automatic}, 
 MaxSteps -> 2 10^4]
(* {0.00158325, 3.17326*10^-6, -0.00026645, -0.00130408} *) 

enter image description here

  • $\begingroup$ Thanks a lot for such a detailed answer! I am getting the first plot but I am not sure why I am still getting errors and warning messages for the Lyapunov Exponents: imgur.com/a/h9YWPag . The code is still running and I am not able to get the main plot. $\endgroup$
    – codebpr
    May 4 at 6:29
  • $\begingroup$ I also got this error for different initial conditions. Do you use the same parameters / ICs as in my answer? In any case, it's because the system is stiff and once it fails, it's a lost cause. You could try adding NDSolveOpts -> {AccuracyGoal -> 16} to see if that helps, or play with other NDSolve options. $\endgroup$
    – Chris K
    May 4 at 16:37
  • $\begingroup$ I have used the same parameters as in your answer. I will try to use other NDSolve options and let you know if I get similar results. $\endgroup$
    – codebpr
    May 4 at 16:49
  • $\begingroup$ OK! My $Version is 13.0.1 for Mac OS X ARM (64-bit) (January 28, 2022), which might account for the different outcome. $\endgroup$
    – Chris K
    May 4 at 16:53
  • 1
    $\begingroup$ The easiest would be to add to LyapunovExponents at the end, by the other ListPlot command, something like Print[ListPlot[Total[Transpose[edat]], Evaluate[Sequence @@ plotopts]]]. $\endgroup$
    – Chris K
    May 6 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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