Return to Answer

References:

Sandri M (1996) Numerical calculation of Lyapunov exponents. The Mathematica Journal 6:78–84. https://library.wolfram.com/infocenter/Articles/2902/

Wolf A, Swift JB, Swinney HL, Vastano JA (1985) Determining Lyapunov exponents from a time series. Physica D: Nonlinear Phenomena 16:285–317. https://doi.org/10.1016/0167-2789(85)90011-9

Update Mar 12, 2022: Fixed mistake in implementation of PlotExponents

LyapunovExponents[eqnsin_List, icsin : ({__Rule} | _Association), nlein_Integer: 0, opts___?OptionQ] := Module[{

(* options *)
tstep, maxsteps, ndsolveopts, logbase, showplot, plotexponents, plotopts,

(* other variables *)
δ, 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[δ[i, j]'[t] == (jac.Table[δ[i, j][t], {i, neq}])[[i]], {j, nle}, {i, neq}]]
];
unks = Join[
  vars,
  Flatten[Table[δ[i, j], {j, nle}, {i, neq}]]
];
ics = Join[
  Table[var[0] == (var /. icsin), {var, vars}],
  Flatten[Table[δ[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[δ[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[δ[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 -> {}};

Update May 3, 2022: Added c, i, j as local variables

Update Mar 12, 2022: Fixed mistake in implementation of PlotExponents

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 *)
δ, 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[δ[i, j]'[t] == (jac.Table[δ[i, j][t], {i, neq}])[[i]], {j, nle}, {i, neq}]]
];
unks = Join[
  vars,
  Flatten[Table[δ[i, j], {j, nle}, {i, neq}]]
];
ics = Join[
  Table[var[0] == (var /. icsin), {var, vars}],
  Flatten[Table[δ[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[δ[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[δ[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 -> {}};

Here's an updated implementation (works with MMA v11) that generalizes Housam Binous & Nasri Zakia's update to Marco Sandri's package and incorporates ideas from @bbgodfrey.

Here's an updated implementation (works with MMA v11) that generalizes Housam Binous & Nasri Zakia's update to Marco Sandri's package and incorporates ideas from @bbgodfrey. The main advantage is that it uses NDSolve rather than Euler's method for solving the differential equations, so it should be more accurate / faster than Sandri's implementation.