Problem with calculation of Lyapunov exponent

Question

The code given by Chris K for Lyapunov Exponent does not work and gives lot of errors for the dynamical system given in equation(2) of this paper

The dynamical equations:

eqns = {x'[t] == 5*(y[t] - x[t]) + 15*y[t]*z[t], 
   y'[t] == -10*(y[t]*Abs[y[t]]) - y[t] + 4*x[t]*z[t], 
   z'[t] == 3.9*z[t] - x[t] y[t]};

with ics:

ics = {x -> 0, y -> 5, z -> 0};

The phase space diagram of the attractor:

The borrowed code for the Lyapunov exponent is given below:

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@NDSolve`ProcessEquations[Flatten[Join[eqns, ics]], unks, t, 
      Evaluate[Sequence @@ ndsolveopts]];
   (*main loop*)
   edat = Table[newstate = First@NDSolve`Reinitialize[state, ics];
     NDSolve`Iterate[newstate, c tstep];
     sol = NDSolve`ProcessSolutions[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 -> 3, PlotOpts -> {}};

LyapunovExponents[eqns, ics, ShowPlot -> True]                                                                                                                                                              

Errors:

How to resolve this type of error and get the desire result?

Please help us out.

Answer 1

In your differential equations, replace Abs with RealAbs.