Lyapunov Exponent of DDE

What is the procedure to compute Lyapunov exponent for Delay Differential Equation using Mathematica code? If we consider the famous Mackey-Glass Equation:

x'[t] == (a x[t - τ]/ (1 + x[t - τ]^c) - b x[t])


where a=0.2, b=0.1, c=10 and τ can be varied from 10 to 40.

• Commented Nov 14, 2015 at 4:09

This page by J.C. Sprott suggests a way to calculate Lyapunov exponents of delay differential equations by approximating them by a large system of ordinary differential equations. Define $x_0(t)=x(t)$ and $x_n(t)=x(t-\tau)$. Then append $n$ equations $dx_i/dt=(n/(2\tau))(x_{i-1}-x_{x+1})$ for $1\leq i\leq n-1$ and $dx_n/dt=(n/\tau)(x_{n-1}-x_n)$. This system can be analyzed using standard methods for calculating Lyapunov exponents of ODEs.

First simulate the original DDE system in the chaotic realm:

a = 0.2; b = 0.1; c = 10; τ = 23;
sol = NDSolve[{x'[t] == a x[t - τ]/(1 + x[t - τ]^c) - b x[t], x[0] == 0.1}, x,
{t, 0, 1000}][[1]];
Plot[x[t] /. sol, {t, 0, 1000}]


Now simulate the approximation by ODEs (here $n=100$):

n = 100;
xp[0] = a x[n][t]/(1 + x[n][t]^c) - b x[0][t];
Do[
xp[i] = n/(2 τ) (x[i - 1][t] - x[i + 1][t]);
, {i, 1, n - 1}];
xp[n] = n/τ (x[n - 1][t] - x[n][t]);

warm = NDSolve[{Table[x[i]'[t] == xp[i], {i, 0, n}],
Table[x[i][0] == 0.1, {i, 0, n}]},Table[x[i], {i, 0, n}], {t, 0, 1000}][[1]];
Plot[x[0][t] /. warm, {t, 0, 1000}]


Looks similar enough, with the discrepancies at the end inevitable given that it's chaotic!

Now define LyapunovExponents as in this answer, then take the final conditions of the warmup as initial conditions for LyapunovExponents. Here I calculate only the largest exponent and take large time steps for each point accumulated, since this is slow enough already (almost 10 minutes).

ics = Table[x[i] -> (x[i][1000] /. warm), {i, 0, n}];
LyapunovExponents[Table[x[i]'[t] == xp[i], {i, 0, n}], ics, 1,
TStep -> 100, ShowPlot -> True]


(* {0.00942851} *)


This is close to the value Sprott got for the same $n=100$ approximation (0.00936), which is itself close to the value he cites as the real answer (0.00956).

• but on the x-axis, it should be delay instead of time. Commented Aug 30, 2018 at 10:41

Following your comment, I have now amended this code to incorporate the parameter estimates you suggested.

  T = 20; c = 10;
y = ParametricNDSolveValue[{Derivative[1][x][t] ==
a*(x[t - T]/(1 + x[t - T]^c)) - x[t]*b,
x[t /; t <= 0] == 0.1}, {x, Derivative[1][x]}, {t, 0, 500}, {a,
b}];

Manipulate[
ParametricPlot[{y[a, b][[1]][t], y[a, b][[2]][t]}, {t, 0, 500},
AxesLabel -> {x, Derivative[1][x]}, AspectRatio -> 1], {{a, 0.1},
0.05, 0.5}, {{b, 0.1}, 0.05, 0.3}]


The x-axis is x[t] and the y-axis s x'[t] and there appears to be oscillations at select range of a and b values.

What is the suitable value for T? The code below may help

  Clear[T, a, b, c]
Manipulate[
Module[{a = 0.2, b = 0.1, c = 10, eqs},
eqs = NDSolve[{Derivative[1][x][t] ==
a*(x[t - T]/(1 + x[t - T]^c)) - x[t]*b,
x[t /; t <= 0] == 0.1}, x,
{t, 0, 500}];
Plot[Evaluate[{x[t], Derivative[1][x][t]}] /. eqs, {t, 0, 500},
PlotRange -> All]], {{T, 15}, 10, 50}]


The Lyapunov exponent is a measure of sensitive dependence on conditions at t=0 and is calculated based on how rapidly two nearby states diverge from each other. So evaluation of this exponent depends on the parameter that you have selected. I just show the evaluation of the Lyapunov Exponent for a simpler case where it exist, such as f(x)=a x(b-x) where a and b are parameters. I set b=1.

lya[a_, b_] :=
Module[{f, ly}, f[x_] := a*x*(b - x);
ly[x_] := Log[Abs[a*(b - 2*x)]]; Mean[ly[NestList[f, 0.05, 100]]]];
Show[Plot[lya[a, 1], {a, 0, 4}, PlotTheme -> "Scientific",
PlotStyle -> Thickness[0.01], PlotRange -> {-3, 1.5}],
ListPlot[
ParallelTable[
Thread[{a, Nest[a*#1*(1 - #1) & , Range[0, 1, 0.01], 500]}], {a,
0, 4, 0.01}]]]


Notice that there are locations where nearby points diverge and the lyapunov exponent swings downwards. This code can be modified to determine the lyapunov exponent for your model (provided you fix all known parameters except the one you wish to test for sensitivity of evolving states). No doubt there are alternative approaches to determine the lyapunov exponent for higher dimensional cases. I have just shown the one-dimensional case here.

• Thank you for your reply. But how to find out Lipunov Exponent? T can be Varied from 10 to 40. and t can be varied from 0 t0 500 Commented Nov 14, 2015 at 4:01
• I have amended the code to better reflect the estimates you suggested. Can you modify the code to evaluate the Lyapunov exponent....for the specific parameter you have in mind? Commented Nov 14, 2015 at 11:03