# System of delay differential equations

I'm trying to solve numerically the following system of delay differential equations

$$\dot{x} = \lambda x(t)+ \omega y(t) - K(t) (x(t)- x(t-\tau)) \\ \dot{y} = -\omega x(t) + \lambda y(t) - K(t) (y(t)- y(t-\tau)) \\ \dot{K} = \gamma ((x(t)-x(t-\tau))(x(t)-2x(t-\tau)+x(t-2\tau)) + (y(t)-y(t-\tau) (y(t)-2y(t-\tau)+y(t-2\tau)))$$

with history functions $x(t)=y(t)=0$ for $t<0$ and $K(t)=0$ for $t\leq2\tau$ and initial conditions $x(0)=0.02$, $y(0)=0$.

I first tried to do it with $x(0)=0$ which should result in $x(t)=y(t)=K(t)=0$ for all $t$ but I already get an error message saying that (at least) one derivative at $t=0$ is non-numerical. Here is my code

lam = 0.5; om = Pi; tau = 1; gam = 1;

NDSolve[{x'[t] == lam x[t] + om y[t] - k[t] (x[t] - x[t - tau]),
y'[t] == -om x[t] + lam y[t] - k[t] (y[t] - y[t - tau]),
k'[t] == gam ((x[t] - x[t - tau]) (x[t] - 2 x[t - tau] +
x[t - 2 tau]) + (y[t] - y[t - tau]) (y[t] - 2 y[t - tau] +
y[t - 2 tau])), x[t /; t <= 0] == 0,
y[t /; t <= 0] == 0, k[t /; t <= 2 tau] == 0}, {x, y, k}, {t, - tau, 40}]

• If you change initial condition k[t /; t <= 2 tau] == 0 to k[t /; t <= 0] == 0 it seems to work. Jun 13, 2018 at 15:58
• Thanks for the suggestion but I would like to keep it that way. Also how can i change the initial value for $x$ to $x(0)=0.02$? Jun 13, 2018 at 18:34
• If you use k[t /; t <= 0] == HeavisideTheta[t - 2 tau]  you get a solution, but all values of y,y,k are zero. If you apply Reduce[{eqs, t == -1}]  there seem to occur condtradictions. Jun 14, 2018 at 7:01

The following works:

First we solve the ordinary differential equation system for $x$ and $y$ up to $t=2\tau$ without considering $k$.

om = Pi; lam = 0.5; tau = 1; gam = 1;
solstart = {x[t], y[t]} /.
NDSolve[{x'[t] ==lam x[t]+om y[t],y'[t]==-om x[t]+lam y[t],
x == 0.5, y == 0}, {x, y}, {t, 0, 2tau}][];


Then we use these solutions as history functions for the delay differential equation system. The whole code:

om = Pi; lam = 0.5; tau = 1; gam = 1;
solstart = {x[t], y[t]} /.
NDSolve[{x'[t]==lam x[t]+om y[t],y'[t]==-om x[t]+lam y[t],
x == 0.5, y == 0}, {x, y}, {t, 0, 2}][];

sol = NDSolve[{x'[t] == lam x[t] + om y[t] - k[t] (x[t] - x[t - tau]),
y'[t] == -om x[t] + lam y[t] - k[t] (y[t] - y[t - tau]),
k'[t] == gam ((x[t] - x[t - tau]) (x[t] - 2 x[t - tau] +
x[t - 2 tau]) + (y[t] - y[t - tau]) (y[t] - 2 y[t - tau] +
y[t - 2 tau])), x[t /; t <= 2tau] == solstart[],
y[t /; t <= 2tau] == solstart[], k[t /; t <= 2tau] == 0},
{x,y,k},{t, 0, 40}][];

Plot[{x[t], y[t]} /. sol // Evaluate, {t, 0, 40}, ImageSize -> Large,
PlotRange -> {{0, 40}, {-1.5, 1.5}}]