3
$\begingroup$

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:

enter image description here

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:

enter image description here

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

Please help us out.

$\endgroup$
1
  • $\begingroup$ Welcome to Mathematica.SE, Arssat! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – Chris K
    Feb 29 at 17:57

1 Answer 1

6
$\begingroup$

In your differential equations, replace Abs with RealAbs.

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Glad it wasn't my fault ;) $\endgroup$
    – Chris K
    Mar 1 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.