Calculating the largest Lyapunov exponent with NDSolve for a damped pendulum

Question

I'm using Mathematica to compute the maximum Lyapunov exponent for a damped pendulum. Ideally, when the damping is zero, the exponent should be zero, and for finite damping, it should be negative. However, even when I set the damping factor $b$ to $0$, I'm obtaining values like $-2.0$. I initially thought that increasing the time step and reducing the initial displacement might enhance the result, but this doesn't seem to be the case. I've also attempted to decrease the maximum step size, but it hasn't proven helpful. I've encountered similar issues when applying this approach to other models.

Can you help me identify where I might be making a mistake? Is my definition even correct?

tmax = 1100;
dellya = 10^(-9);
dt = 0.01;
steps = tmax/dt;
initial[k_, b_] := NDSolve[{phi1''[t] == -k*phi1[t] - b*phi1'[t], phi1[0] == 0.5, 
  phi1'[0] == 0.0}, {phi1[t]}, {t, 0, tmax}, 
  Method -> "ExplicitRungeKutta", MaxStepSize -> 0.1, MaxSteps -> 10^14];
displaced[k_, b_] := NDSolve[{phi1''[t] == -k*phi1[t] - b*phi1'[t], 
  phi1[0] == 0.5 - dellya, phi1'[0] == 0.0}, {phi1[t]}, {t, 0, 
  tmax}, Method -> "ExplicitRungeKutta", MaxStepSize -> 0.1, MaxSteps -> 10^14];
sol1 = initial[2.0, 0.0];
sol2 = displaced[2.0, 0.0];
phi1sol[t_] := Evaluate[phi1[t] /. sol1];
phi2sol[t_] := Evaluate[phi1[t] /. sol2];
phiplot1 = Flatten@Table[phi1sol[t], {t, 0, tmax, dt}];
phiplot2 = Flatten@Table[phi2sol[t], {t, 0, tmax, dt}];
doval = Sqrt[(phiplot1[[1]] - phiplot2[[1]])^2];
lypunovcdti = Table[Sqrt[(phiplot1[[i]] - phiplot2[[i]])^2], {i, 1, Length@phiplot2}];
maximallyapunov = (1/tmax)*Total[Table[Log[(lypunovcdti[[i]]/doval)], 
  {i, 1, Length[lypunovcdti]}]];

codebpr · Answer 1 · 2024-01-24T11:20:26.547

You can use the algorithm designed by @Chris K for obtaining the lyapunov exponents:

GramSchmidt[w_?MatrixQ] := 
 Module[{v = ConstantArray[0, Length[w]]}, 
  Table[v[[n]] = 
    w[[n]] - 
     Sum[(v[[i]] . w[[n]]/v[[i]] . v[[i]])*v[[i]], {i, n - 1}], {n, 
    Length[w]}];
  v]
LyapunovExponents[eqnsin_List, icsin : ({__Rule} | _Association), 
   nlein_Integer : 0, opts___?OptionQ] := 
  Module[{(options)tstep, maxsteps, ndsolveopts, logbase, showplot, 
    plotexponents, plotopts,(iterators)c, i, 
    j,(other variables)[Delta], neq, nle, vars, rhs, jac, eqns, 
    unks, ics, cum, res, edat, state, newstate, sol, W, 
    norms},(parse options)
   tstep = Evaluate[
     TStep /. Flatten[{opts, Options[LyapunovExponents]}]];
   maxsteps = 
    Evaluate[
     MaxSteps /. Flatten[{opts, Options[LyapunovExponents]}]];
   ndsolveopts = 
    Evaluate[
     NDSolveOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
   logbase = 
    Evaluate[LogBase /. Flatten[{opts, Options[LyapunovExponents]}]];
   showplot = 
    Evaluate[
     ShowPlot /. Flatten[{opts, Options[LyapunovExponents]}]];
   plotexponents = 
    Evaluate[
     PlotExponents /. Flatten[{opts, Options[LyapunovExponents]}]];
   plotopts = 
    Evaluate[
     PlotOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
   neq = Length[eqnsin];
   If[nlein == 0, nle = neq, 
    nle = nlein];(how many exponents)(extract vars and right hand 

sides from eqnsin)vars = eqnsin[[All, 1, 0, 1]];
   rhs = eqnsin[[All, 2]];
   (jacobian matrix)jac = D[rhs, {Replace[vars, {x_ -> x[t]}, 1]}];
   eqns = 
    Join[eqnsin, 
     Flatten[Table[[Delta][i, j]'[
         t] == (jac . Table[[Delta][i, j][t], {i, neq}])[[i]], {j, 
        nle}, {i, neq}]]];
   unks = 
    Join[vars, Flatten[Table[[Delta][i, j], {j, nle}, {i, neq}]]];
   ics = 
    Join[Table[var[0] == (var /. icsin), {var, vars}], 
     Flatten[Table[[Delta][i, j][0] == 
        IdentityMatrix[neq][[i, j]], {j, nle}, {i, neq}]]];
   cum = Table[0, {nle}];
   state = 
    First@NDSolveProcessEquations[Flatten[Join[eqns, ics]], unks, t, Evaluate[Sequence @@ ndsolveopts]]; (*main loop*) edat = Table[newstate = First@NDSolveReinitialize[state, ics];
     NDSolveIterate[newstate, c tstep]; sol = NDSolveProcessSolutions[newstate];
     W = 
      GramSchmidt[
       Evaluate[
        Table[[Delta][i, j][c  tstep], {j, nle}, {i, neq}] /. sol]];
     norms = Map[Norm, W];
     (update running vector magnitudes)
     cum = cum + Log[logbase, norms];
     ics = 
      Join[Table[
        var[c  tstep] == (var[c  tstep] /. sol), {var, vars}], 
       Flatten[Table[[Delta][i, j][
           c  tstep] == (W/norms)[[j, i]], {j, nle}, {i, neq}]]];
     cum/(c  tstep), {c, maxsteps}];
   If[showplot, 
    Print[ListPlot[Transpose[edat][[1 ;; plotexponents]], 
      Evaluate[Sequence @@ plotopts]]]];
   Return[cum/(maxsteps  tstep)]];
Options[LyapunovExponents] = {NDSolveOpts -> {}, TStep -> 1, 
   MaxSteps -> 10^4, LogBase -> E, ShowPlot -> False, 
   PlotExponents -> 1, PlotOpts -> {}};

Now coming to your equation, first we need to convert it into two first order differential equations so that we can apply the above algorithm on it:

originaleqn = {phi1''[t] == -k  phi1[t] - b  phi1'[t]};
finaleqns = 
 Flatten@Join[
   Solve[originaleqn, 
      phi1''[t]] /. {phi1''[t] -> Derivative[1][dphi1][t], 
      Derivative[1][phi1][t] -> dphi1[t]} /. (lhs_ -> rhs_) -> (lhs ==
        rhs), {phi1'[t] == dphi1[t]}]

We then state your constant values and use your initial conditions to solve:

k = 2.0; b = 0.0; tmax = 1100;
sol = NDSolve[{finaleqns, {phi1[0] == 0.5, dphi1[0] == 0.0}}, {phi1, 
     dphi1}, {t, 0, tmax}, MaxStepSize -> 0.01][[1]];

We obtain the following Parametric plot from the solution:

ParametricPlot[{phi1[t], dphi1[t]} /. sol, {t, 900, tmax}, 
 PlotRange -> All, AspectRatio -> 1]

Now we use the final values to compute the lyapunov exponent:

ics = {phi1 -> (phi1[tmax] /. sol), dphi1 -> (dphi1[tmax] /. sol)};
LyapunovExponents[finaleqns, ics, ShowPlot -> True, 
 PlotExponents -> 2, 
 PlotOpts -> {FrameLabel -> {"Steps", "[Lambda]"}, 
   PlotRange -> {-0.00004, 0.00004}, PlotTheme -> "Detailed", 
   ImageSize -> 800, LabelStyle -> Directive[Black, 20]}, 
 MaxSteps -> 10^5]

The value of two lyapunov exponents that we obtained from above procedure is:

(* {-8.69245*10^-7, 7.5945*10^-7} *)

which are both approximately equal to $0$.

Results for finite damping

These are the results that I obtained for a finite damping b=0.5:

Special thanks to @Chris K for developing this useful algorithm! — codebpr, Jan 24 '24 at 08:37

Calculating the largest Lyapunov exponent with NDSolve for a damped pendulum

1 Answers1