3
$\begingroup$

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]}]
Improve this question
$\endgroup$
4
  • $\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
    Commented May 2, 2022 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
    Commented May 3, 2022 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
    Commented May 3, 2022 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
    Commented May 4, 2022 at 6:25

1 Answer 1

Reset to default
4
$\begingroup$

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

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

Improve this answer
$\endgroup$
13
  • $\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
    Commented May 4, 2022 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
    Commented May 4, 2022 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
    Commented May 4, 2022 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
    Commented May 4, 2022 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
    Commented May 6, 2022 at 16:55

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.