I'm struggling now for several weeks with a very bizarre problem with a system of delay differential equations.

First, here the system:

$$\dot a = 1 - \Theta(b(t-\tau)-\kappa) \,- a(t) \\ \dot b = \,\,\,\,\,\,\,\,\,\,\,\Theta(a(t-\tau)-\kappa) - b(t),$$

where $$\Theta: \Bbb R \rightarrow \{ 0,1 \},\,\,x \mapsto \begin{cases} 1, x \ge 0 \\ 0, x < 0 \end{cases}$$ is the Heaviside function, $\tau$ a delay and $0 \le \kappa \le 1$ a constant.

And here a visualisation of the numerically calculated periodic dynamics of $a$ and $b$ (please ignore the variable $c$):

enter image description here

The horizontal dashed line marks the value of $\kappa$, here $\kappa \approx 0.364$ was used. Furthermore, $\tau = 0.89$ was used as well as the initial conditions $a(0) = 0$ and $b(0) = 0$ (the plot shows the dynamics after the transient oscillations).

We want to solve the system of delay differential equations given above analytically by parting the interval $[0, T]$ into appropriate sections.

First, let's have a look at the dynamics of $a$:

  • $0 \leq t \leq t_a^{max}$: \begin{align} \partial_t a_1\left(t\right) &= 1 - a_1\left(t\right) \nonumber \\ \Rightarrow a_1\left(t\right) &= \left(a_{min} - 1\right) \mathrm{e}^{-t} + 1 \end{align}

  • $t_a^{max}\leq t \leq T$: \begin{align} \partial_t a_2\left(t\right) &= - a_2\left(t\right) \nonumber \\ \Rightarrow a_2\left(t\right) &= a_{max} \mathrm{e}^{-\left(t-t_a^{max}\right)} \end{align}

The continuity conditions are given by \begin{align} a_1 \left(0\right) &= a_2\left(T\right) \nonumber \\ a_1\left(t_a^{max}\right) &= a_2\left(t_a^{max}\right). \end{align}

We get \begin{align} a_{min} &= \frac{\mathrm{e}^{t_a^{max}-T} - \mathrm{e}^{-T}}{1-\mathrm{e}^{-T}} \nonumber \\ a_{max} &= \frac{1- \mathrm{e}^{-t_a^{max}}}{1-\mathrm{e}^{-T}}. \end{align}

Now, let's have a look at the dynamics of $b$:

Analogous to the dynamics of $a$, we get

  • $0 \le t \le t_b^{min}$: \begin{align} b_0(t) = b_r e^{-t} \end{align}

  • $t_b^{min} \le t \le t_b^{max}$: \begin{align} b_1\left(t\right) = \left(b_{min} - 1\right) \mathrm{e}^{-\left(t-t_b^{min}\right)} + 1 \nonumber \end{align}

  • $t_b^{max} \le t \le T$: \begin{align} b_2\left(t\right) = b_{max} \mathrm{e}^{-\left(t-t_b^{max}\right)} \end{align}

and with the continuity conditions \begin{align} b_0(0) &= b_2(T) \\ b_1 \left(t_b^{min}\right) &= b_2\left(t_b^{min}\right) \nonumber \\ b_1\left(t_b^{max}\right) &= b_2\left(t_b^{max}\right) \end{align} we conclude \begin{align} b_{min} &= \frac{\mathrm{e}^{t_b^{max} - t_b^{min} -T} - \mathrm{e}^{-T}}{1-\mathrm{e}^{-T}} \nonumber \\ b_{max} &= \frac{1- \mathrm{e}^{t_b^{min}-t_b^{max}}}{1-\mathrm{e}^{-T}}. \end{align}

If we calculate now the times $t_a^{max}$ and $T$ for a given delay $\tau$, we see something very bizarre (the calculation is not shown here because it's very long, but we are very confident, really very confident, that there was no error made in this calculation). We get the following dynamics (a in blue, b is not shown):

enter image description here

Obviously, the curve of $a(t)$ is horizontally mirrored in compare to the curve of $a(t)$ shown in the image above which we get from the numerical solution.

How is it possible that we now get this curve instead of a similar curve to the one in the first image? Did we make any wrong assumptions? Can anybody solve this mystery?

Some other hints that may help to solve the mystery:

In the analytical solution it holds $a_{max} < a_{min}$ and $b_{min} > b_{max}$. This comes along with the fact that the times $t_a^{max}$ and $T$ are negative (there is no positive solution for them).

  • $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    May 3, 2017 at 18:29
  • $\begingroup$ Well, I think that is a good idea. Is it allowed to recreate this question there or do I have to delete it here before? $\endgroup$
    – Peter123
    May 3, 2017 at 19:18
  • 1
    $\begingroup$ It's unclear to me: is it mirrored, or shifted to the left? The two graphs' time axes don't match, so it's hard to see. $\endgroup$
    – Kirill
    May 3, 2017 at 20:36
  • $\begingroup$ It's really mirrored. You can see that also by the yellow curve (that's the solution for another variable which is coupled to a). $\endgroup$
    – Peter123
    May 3, 2017 at 21:06
  • $\begingroup$ If the plot shows the transient dynamics, how come the red and blue curve not starting at 0? $\endgroup$
    – Wrzlprmft
    May 5, 2017 at 17:05

1 Answer 1


First of all, you have no exact mirroring, which can be clearly seen when superimposing your plots:

superimposed plots

I suggest to consider the following to get a better idea of the possible mistake:

  • Consider different sets of parameters, in particular ones with a less symmetric dynamics, e.g., $κ=0.2$ and $τ=0.1$. If there is a disagreement between numerical and analytical results, it should be much clearer then.

  • I do not know what integrator you used, but if it relies on interpolating and adapts step sizes, the discontinuity of the Heaviside function may cause problems. I had to replace it with a steep sigmoidal function to numerically integrate the system with adaptive step sizes. (Whether my results agree with yours, I cannot say as you do not give any initial conditions or $τ$.)

  • $\begingroup$ Thank you very much for your answer. I'll check that. I also used a steep sigmoidal function: $f(x) = \frac{x^n}{\kappa^n+x^n}$. I used $\tau = 0.89$ and $a(0) = 0$ as well as $b(0) = 0$ as initial conditions. $\endgroup$
    – Peter123
    May 5, 2017 at 12:07
  • $\begingroup$ Sure, sorry. Was a typo. $\endgroup$
    – Peter123
    May 5, 2017 at 12:13
  • $\begingroup$ @Peter123: You might want to edit that information into your question. $\endgroup$
    – Wrzlprmft
    May 5, 2017 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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