Mysterious Mirroring in Analytical Solution of a delay differential equation (DDE)

Question

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$):

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):

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).

Answer 1

First of all, you have no exact mirroring, which can be clearly seen when superimposing your 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 $τ$.)