# Lyapunov exponent of Delay Differential Equation

I have a Delay Differential Equation as given below

b = 1;
V0 = 300;
I0 = 0.001;
N1 = 7/4;
p1 = (π*b*N1)/V0;
a1 = 1.001; a2 = 0.123; a3 = -3.622*10^-3; b1 = 0.001959; b2 = 0.031; \
b3 = 0.003241; G = 0.5*10^-5; Ω = 0; C1 = p1/2000;
f = 1*10^3;
Is = 0.001;
τ = 10;
NL = 1000;

sol1[t_] =
NDSolve[{x'[t] - (I0*p1)/(
2*C1)*((a1*x[t - τ] - G/(p1*I0)*x[t]) -
3/4*x[t - τ]^3*a2 - 5/8*a3*x[t - τ]^5) -
Is*p1*Cos[y[t]] == 0,
y'[t] - Ω - (I0*p1)/(
2*C1*x[t])*(b1*x[t - τ] + 3/4*x[t - τ]^3*b2 -
5/8*b3*x[t - τ]^5) + (Is*p1)/(2*C1*x[t])*Sin[y[t]] == 0,
x[0] == 0.003, y[0] == 0.001}, {x, y}, {t, 0, NL}]

I want to compute the Lyapunov exponent for this system by varying Tau from 3 to 22. For This I have used the following program.

b = 1;
V0 = 300;
I0 = 0.001;
N1 = 7/4;
p1 = (π*b*N1)/V0;
a1 = 1.001; a2 = 0.123; a3 = -3.622*10^-3; b1 = 0.001959; b2 = 0.031; \
b3 = 0.003241; G = 0.5*10^-5; Ω = 0; C1 = p1/2000;
f = 1*10^3;
Is = 0.001;
τ = 10;
NL = 1;

deq1 = (I0*p1)/(
2*C1)*((a1*x1[t - τ] - G/(p1*I0)*x1[t]) -
3/4*x1[t - τ]^3*a2 - 5/8*a3*x1[t - τ]^5) +
Is*p1*Cos[y1[t]];
deq2 = Ω + (I0*p1)/(
2*C1*x1[t])*(b1*x1[t - τ] + 3/4*x1[t - τ]^3*b2 -
5/8*b3*x1[t - τ]^5) - (Is*p1)/(2*C1*x1[t])*
Sin[y1[t]];
x10 = 0.003;
y10 = 0.001;
dx0 =  10^-8;
tin = 0;
tfin = 201;
tstep = NL;
acc = 12;
lcedata = {};
sum = 0;
d0 = Sqrt[(x10)^2 + (y10)^2];
For[i = 1, i < tfin/tstep, i++,
sdeq = {x1'[t] == deq1, y1'[t] == deq2, x1[0] == x10, y1[0] == y10};
sol = NDSolve[sdeq, {x1[t], y1[t]}, {t, 0, tstep},
MaxSteps -> Infinity, Method -> "StiffnessSwitching",
PrecisionGoal -> acc, AccuracyGoal -> acc];

xx1[t_] = x1[t] /. sol[[1]];
yy1[t_] = y1[t] /. sol[[1]];
d1 = Sqrt[(xx1[tstep])^2 + (yy1[tstep])^2];
sum += Log10[d1/d0];
dlce = sum/(tstep*i);
AppendTo[lcedata, {tstep*i, Log10[dlce]}];
w1 = (xx1[tstep])*(d0/d1);
w2 = (yy1[tstep])*(d0/d1);
x10 = xx1[tstep];
y10 = yy1[tstep];
x20 = x10 + w1;
y20 = y10 + w2;
i = i++;
If[Mod[tstep*i, 1] == 0,
Print[" For t = ", tstep*i, " , ", " LCE = ", Log10[dlce]]]]
S0 = ListLinePlot[{lcedata}, PlotRange -> Automatic,
AxesLabel -> {"t", "log10(LCE)"},
AxesStyle -> Directive["Black", 13], GridLines -> Automatic]

This program gives me some values of the Lyapunov exponent for a given Tau (say Tau=10), but if I change the value of Tau the Lyapunov exponent values are not changing. Please help me to modify the program so that one can get different Lyapunov exponent values for different Tau.

• I'm getting after evaluating sol1[t_]: StringForm::sfr: "Item 2 requested in \""Delayed time 1 = 2 computed at 3 = 4 did" not evaluate to a real number. Sep 5, 2016 at 16:48
• Tau should be \[Tau] in the initialization (probably a typo in the question). More importantly, NDSolve initial conditions are not properly configured for a delayed ODE, which causes the result to be independent of \[Tau]. Sep 5, 2016 at 19:47
• I'm voting to close this question as off-topic because it's too localized and unlikely to help future visitors. Jan 4, 2017 at 11:23
• For delay equations, not only the variables but also all the variable values at previous time delay history (if delay is 1 and integration step 0.1, there shall be 10 values of variable at previous time steps) need to be considered as system variables for calculating the growth vector. Also in the next iteration loop, initial condition as well as delay history has to be normalized. Finally, the system equations need to be linearized (Jacobian) before integrating inside a loop for better results. Sep 6, 2017 at 10:32

Let's pick one $$\tau$$ value, say $$\tau=18$$. First, simulate the DDEs:

τ = 18;
sol = NDSolve[{x'[t] == (I0*p1)/(2*C1)*
((a1*x[t - τ] - G/(p1*I0)*x[t]) - 3/4*x[t - τ]^3*a2 - 5/8*a3*x[t - τ]^5)
+ Is*p1*Cos[y[t]],
y'[t] == Ω + (I0*p1)/(2*C1*x[t])*(b1*x[t - τ] + 3/4*x[t - τ]^3*b2 -
5/8*b3*x[t - τ]^5) - (Is*p1)/(2*C1*x[t])*Sin[y[t]],
x[0] == 0.003, y[0] == 0.001}, {x, y}, {t, 0, NL}];
Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, NL}]

I almost gave up at this point, because it doesn't look like any attractor is reached (particularly for y[t]). Then I noticed that y[t] only shows up as an argument of Sin[y[t]] and Cos[y[t]], which are periodic, so I guess this might be considered as a circular phase space in y[t].

Plot[Evaluate[x[t] /. sol], {t, NL - 200, NL}]

The dynamics of x[t] look decent. Let's proceed (at our own peril?)

Following J.C. Sprott's suggestion and the approach I used in this answer, we can replace the DDEs with a large set of ODEs and then calculate its Lyapunov exponents.

n = 40;
xp[0] = (I0*p1)/(2*C1)*((a1*x[n][t] - G/(p1*I0)*x[0][t]) -
3/4*x[n][t]^3*a2 - 5/8*a3*x[n][t]^5) + Is*p1*Cos[y[0][t]];
yp[0] = Ω + (I0*p1)/(2*C1*x[0][t])*(b1*x[n][t] +
3/4*x[n][t]^3*b2 - 5/8*b3*x[n][t]^5) - (Is*p1)/(2*C1*x[0][t])* Sin[y[0][t]];
Do[
xp[i] = n/(2 τ) (x[i - 1][t] - x[i + 1][t]);
yp[i] = n/(2 τ) (y[i - 1][t] - y[i + 1][t])
, {i, 1, n - 1}];
xp[n] = n/τ (x[n - 1][t] - x[n][t]);
yp[n] = n/τ (y[n - 1][t] - y[n][t]);

Warm up to get on attractor:

warm = NDSolve[Flatten[Join[
Table[{x[i]'[t] == xp[i], y[i]'[t] == yp[i]}, {i, 0, n}],
Table[{x[i][0] == 0.003, y[i][0] == 0.001}, {i, 0, n}]
]], Flatten[Table[{x[i], y[i]}, {i, 0, n}]], {t, 0, 1000}][[1]];

Plot[Evaluate[x[0][t] /. warm], {t, 0, 1000}]

Use final conditions as initial conditions for LyapunovExponents function from this answer. We'll get only the largest exponent to save time.

ics = Flatten[Table[{
x[i] -> (x[i][1000] /. warm),
y[i] -> (y[i][1000] /. warm)}, {i, 0, n}]];

LyapunovExponents[Flatten[Table[{x[i]'[t] == xp[i], y[i]'[t] == yp[i]}, {i, 0, n}]],
ics, 1, ShowPlot -> True]

(* {0.00857852} *)

Looping this across $$\tau$$'s is left as an exercise for the reader :)

EDIT 4/29/2020: In the comments, Udichi asked about a second-order system that had been converted to two first-order equations. This illustrates that only variables with delays need to be converted to a system.

ϵ = 0.2; τ = 9; ϕ = -(π/4); β = 3;
sol = NDSolve[{
y'[t] == x[t],
x'[t] == -x[t] - ϵ y[t] + β Cos[x[t - τ] + ϕ]^2,
y[0] == 0.1, x[0] == 0.001}, {y, x}, {t, 0, 1000}][[1]];
Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 1000}, PlotPoints -> 100]

Looks chaotic to the naked eye. Converting x[t] to n = 40 equations at lags up to τ and applying LyapunovExponents to calculate the largest exponent:

n = 40;
xp[0] = -x[0][t] - ϵ y[t] + β Cos[x[n][t] + ϕ]^2;
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]);
yp = x[0][t];

ics = Flatten[Join[
Table[x[i] -> 0.001, {i, 0, n}],
{y -> 0.1}
]];

(* 44 sec *)
LyapunovExponents[
Flatten[Join[Table[x[i]'[t] == xp[i], {i, 0, n}], {y'[t] == yp}]],
ics, 1, ShowPlot -> True, TStep -> 100, MaxSteps -> 1000]

{0.0048002}

• I have a doubt, If in the above-coupled equation, x'[t]==y[t] and y'[t] is same as described above, then what to do with xp[i] and xp[n]? I mean the equation is coupled but only y'[t] is DDE ,not x'[t]. Then how to deal with the problem? could you please give some insight? Apr 26, 2020 at 19:13
• @Udichi I think you'd need the same approach, although when I substituted x'[t]==y[t] it did not produce a nice attractor. But your comment made me realize that my answer was inefficient: since only y[t] is used, there is no need for all the delayed versions, y[i][t]. If that's useful for your to see, let me know and I'll edit my answer. Apr 28, 2020 at 1:44
• thanks for your response. Your script for the above equation is nice, but in many problems, we have 2nd order DDE. We can decompose it into two equations, out of which one may be a DDE and the other may not be. what I tried, xp [i]= y[i] [t] and xp[n] = y[n] [t] , but I did not get success.Is this approach correct? Stay safe. Apr 28, 2020 at 4:43
• @Udichi If you add an example of a second order DDE to your question, we could take a look. Apr 28, 2020 at 14:46
• @Udichi Check out my edit ^^ Apr 30, 2020 at 2:01

The statement Tau = 10; should be τ = 10; (probably a typo in the question). With this change made, NDSolve complains about initial conditions. Use instead,

sol1 = Flatten @ NDSolve[{x'[t] - (I0*p1)/(2*C1)*((a1*x[t - τ] - G/(p1*I0)*x[t]) -
3/4*x[t - τ]^3*a2 - 5/8*a3*x[t - τ]^5) - Is*p1*Cos[y[t]] == 0,
y'[t] - Ω - (I0*p1)/(2*C1*x[t])*(b1*x[t - τ] + 3/4*x[t - τ]^3*b2 -
5/8*b3*x[t - τ]^5) + (Is*p1)/(2*C1*x[t])*Sin[y[t]] == 0,
x[t /; t <= 0] == 0.003, y[t /; t <= 0] == 0.001}, {x, y}, {t, 0, NL}]

Results for τ = 10 are

ParametricPlot[{x[t], y[t]} /. sol1, {t, 0, NL}, AspectRatio -> 1]

and for τ = 2

Thus, the ODE solver now is working properly. Unfortunately, the second half of the code has numerous errors and produces no results. I do not have time to debug this now.

• I have no problem with the parametric Plot. Only I need to find the LE of this differential equation. I have modified the question Sep 6, 2016 at 11:07