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.
$\begingroup$
$\endgroup$
-
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
$\begingroup$
$\endgroup$
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]
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]
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, withnan integer. - The number of steps that must be saved for
Nestto process isCeiling[n(tdelay-dt/2)]. - The offset,
dt/2arises 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])
use
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]
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
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.
