# Convergence plot of Lyapunov exponents

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],
Derivative[1][OverTilde[c][0]][t],
Derivative[1][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]}]

• 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? May 2 at 20:06
• @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) May 3 at 8:06
• Can you get even NDSolve to work on the model? Before calculating the Lyapunov exponents, you should be able to solve the dynamics. May 3 at 15:29
• Actually, I used the above equations to plot a Poincare section using NDSolve and WhenEvent, hence I assumed it to be working. May 4 at 6:25

## 1 Answer

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[
Solve[
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)};


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} *)


• 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. May 4 at 6:29
• 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. May 4 at 16:37
• 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. May 4 at 16:49
• OK! My \$Version is 13.0.1 for Mac OS X ARM (64-bit) (January 28, 2022), which might account for the different outcome. May 4 at 16:53
• 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]]]. May 6 at 16:55