I would like to know if it is possible in Mathematica to simulate a differential equation with varying delay for example the logistic equation: $$x'(t)=x(t)(1-x(t-\tau(t)))$$ where for example $\tau(t)=\sin^2(t)$. It is mentioned that the NDSolve command only allows constant delays.

  • 2
    $\begingroup$ what is the initial condition? $\endgroup$ – Nasser Jul 3 '17 at 23:09
  • 1
    $\begingroup$ I am not aware of any Mathematica function that can solve this ODE, but one probably could be written without too much difficulty. $\endgroup$ – bbgodfrey Jul 4 '17 at 4:17
  • $\begingroup$ Could you transform it into an equation with fixed delay somehow? $\endgroup$ – Chris K Jul 4 '17 at 7:02

As noted in my earlier comment, I am unaware of an existing Mathematica function that can solve the variable delay ODE in the question. Certainly, NDSolve objects, when asked to solve it. So, here is a somewhat rudimentary function written to solve the ODE. But, it does work, and quickly. To begin, here is a problem that NDSolve can handle.

s = x /. Flatten@NDSolve[{x'[t] == x[t] (1 - x[t - 1.05]), 
    x[t /; t <= 0] == Cos[t]}, x, {t, 0, 20}];
pnds = Plot[s[t], {t, 0, 20}, PlotRange -> All]

enter image description here

We now construct a simple function to reproduce this result. Suppose that we have an ODE,

x'[t] == c x[t]

Then, a finite difference approximation to it is

Flatten@Solve[(xp - xm)/dt == c (xp + xm)/2, xp] // Simplify
(* {xp -> -(((2 + c dt) xm)/(-2 + c dt))} *)

Use this as the basis for a solution employing Nest with time-step 1/10.

xx = N@Table[Cos[t], {t, -1, 0, 1/10}];
Last /@NestList[(xold = First[#]; Append[Rest[#], Last[#] (2 + (1 - xold)/10)/
    (2 - (1 - xold)/10)]) &, xx, 200];
pnst = ListPlot[%, DataRange -> {0, 20}, PlotRange -> All]

enter image description here

Superimposing these two curves with Show[pnst, pnds] indicates that they are identical to the eye. Before proceeding, we offer a few observations.

  • The time step must be of the form 1/n, with n an integer.
  • The number of steps that must be saved for Nest to process is Ceiling[n(tdelay-dt/2)].
  • The offset, dt/2 arises from the fact that the finite difference template above is centered at a half-integer time step.

For a delay that varies with time, t must be specified in the calculation. This can be done by using Fold[..., Table[t, {t, 0, 20, 1/10}] instead of Nest[..., 200]. For instance, to solve an only slightly variable delay,

x'[t] == x[t] (1 - x[t - 1.05 + Sin[t]^2/10])


Last /@ FoldList[(xold = Interpolation[#1, 1 + 1 (1 - Sin[#2]^2)]; 
    Append[Rest[#1], Last[#1] (2 + (1 - xold)/10)/(2 - (1 - xold)/10)]) &, xx, 
    Table[t, {t, 0, 20, 1/10}]];
pnst = ListPlot[%, DataRange -> {0, 20}, PlotRange -> All]

enter image description here

which differs only slightly from the second plot in this answer, as expected. For a highly variable delay, for instance,

x'[t] == x[t] (1 - x[t - 1.05 + 9 Sin[t]^2/10])

replace Interpolation[#1, 1 + 1 (1 - Sin[#2]^2)] by Interpolation[#1, 1 + 9 (1 - Sin[#2]^2)] in the previous code block to obtain

enter image description here

in which the oscillation damps much more rapidly. By way of a warning, the case Interpolation[#1, 1 + 10 (1 - Sin[#2]^2)] produces a curve that is not quite accurate at small t, because xold is too close to xm in the original template. Incidentally, generalization to values of dt not equal to 1/n should not be difficult.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.