# 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. – corey979 Sep 5 '16 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]. – bbgodfrey Sep 5 '16 at 19:47
• I'm voting to close this question as off-topic because it's too localized and unlikely to help future visitors. – xzczd Jan 4 '17 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. – mukund Sep 6 '17 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, TMax -> 10^2, ShowPlot -> True]

(* {0.00857852} *)

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

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 – Udichi Sep 6 '16 at 11:07