Bifurcation of the dde $x'(t) = r x(t-d) - h x(t)$ wich describes a population model with delay

Question

Consider the following population model: $x'(t) = r x(t-d) - h x(t)$, with growth rate $r > 0$ and death rate $h > 0$, where only adults that have reached the age $d > 0$ can reproduce.

I have to draw a bifurcation diagram. For $d = 0$ it looks easy but I don't know how to do it for $d > 0$ since it's my first encounter with a DDE. Can someone help me?

My approach:

The first thing to do would be to rescale the time so the delay becomes 1 and have all parameters outside $x$. we can take $a = \frac{t}{d}$, it'll probably look something like this: $\frac{1}{d}x′(a) = ... - h y(a)$, where $y(a) = x(t)$. I don't know how the term $r x(t-d)$ should look. I can also assume we seek solutions of the form $ce^{\lambda t}$, so I can derive the characteristic function.

Answer 1

If I understand it correctly, your problem is solving the DDE. To do this you can simply substitute $x\left( t \right) = e^{\lambda \cdot t}$ (aka $x'\left( t \right) = \lambda \cdot e^{\lambda \cdot t}$) - reducing the DDE to it's characteristic polynomial - and solve for $\lambda$ via the Lambert-W Function: \begin{align*} x'\left( t \right) &= r \cdot x\left( t - d \right) - h \cdot x\left( t \right) \tag{$x\left( t \right) := e^{\lambda \cdot t}$}\\ \lambda \cdot e^{\lambda \cdot t} &= r \cdot e^{\lambda \cdot \left( t - d \right)} - h \cdot e^{\lambda \cdot t}\\ \lambda \cdot e^{\lambda \cdot t} &= r \cdot e^{\lambda \cdot t} \cdot e^{-\lambda \cdot d} - h \cdot e^{\lambda \cdot t} \tag{$\div e^{\lambda \cdot t}$}\\ \lambda &= r \cdot e^{-\lambda \cdot d} - h \tag{$+ h$}\\ \lambda + h &= r \cdot e^{-\lambda \cdot d} \tag{$\cdot \left( d \cdot e^{d \cdot \lambda} \cdot e^{d \cdot h} \right)$}\\ \left( d \cdot \lambda + d \cdot h \right) \cdot e^{d \cdot \lambda + d \cdot h} &= d \cdot r \cdot e^{d \cdot h} \tag{$\operatorname{W}_{k}\left( \cdot \right)$}\\ d \cdot \lambda + d \cdot h &= \operatorname{W}_{k}\left( d \cdot r \cdot e^{d \cdot h} \right) \tag{$- \left( d \cdot h \right)$}\\ d \cdot \lambda &= \operatorname{W}_{k}\left( d \cdot r \cdot e^{d \cdot h} \right) - d \cdot h\tag{$\div d$}\\ \lambda &= \frac{\operatorname{W}_{k}\left( d \cdot r \cdot e^{d \cdot h} \right) - d \cdot h}{d}\\ \end{align*} where $k \in \mathbb{Z}$ and $\operatorname{W}_{k}$ is the Lambert-W Function.

So ($\text{c}_{k}$ and $\text{c}$ are arbitrary complex constants): $$\fbox{$ \begin{align*} x\left( t \right) &= \text{c} \cdot e^{\frac{\operatorname{W}_{k}\left( d \cdot r \cdot e^{d \cdot h} \right) - d \cdot h}{d} \cdot t}\\ \end{align*} $} \tag{1}$$ or $$\fbox{$ \begin{align*} x\left( t \right) &= \sum\limits_{k = -\infty}^{\infty}\left[ \text{c}_{k} \cdot e^{\frac{\operatorname{W}_{k}\left( d \cdot r \cdot e^{d \cdot h} \right) - d \cdot h}{d} \cdot t} \right]\\ \end{align*} $} \tag{2}$$ because there are infinitely many solutions and they are all homogeneous.

Of course, for your scenario not every $k$ and $c$ makes sense, because a population number that is a complex non-real number makes not that much sense. $\text{c}$ could be interpreted as the initial stock / amount of people and the most realistic choice for $k$ would probably be $k = 0$ (often referred to as the main branch) or $\text{c}_{k} = 0 \ne k \wedge \text{c}_{0} := \text{c}$.

---

That would be a plot for $x\left( t \right)$ where $\text{c} := 100$, $r := 1.5$, $d := 20$ and $h = 0.1$ ($t \in \left[ -1,\, 10 \right]$):

That would be a plot with the parameters for $x_{r}\left( t \right)$ where $r$ goes from $0$ to $2$ in $0.2$ steps and $t \in \left[ -1,\, 20 \right]$:

Mathematica-Code:

c := 100; d := 20; h := 0.1;
Subscript[x, r_][t_] := c * Exp[((LambertW[0, d * r * Exp[d * h]] - d * h) / d) * t]
Plot[Subscript[x, 1.5][t], {t, -1, 10}]

and

Plot[{Subscript[x, 0][t], Subscript[x, 0.2][t], Subscript[x, 0.4][t], Subscript[x, 0.6][t], Subscript[x, 0.8][t], Subscript[x, 1][t], Subscript[x, 1,2][t], Subscript[x, 1.4][t], Subscript[x, 1.6][t], Subscript[x, 1.8][t], Subscript[x, 2][t]}, {t, -1, 20}, PlotLegends -> Table[Subscript["x", N[r / 5]], {r, 0, 10}]]

---

You should be able to work with that. You could also derive solutions on other ways like seeing that there is the trivial solution $x\left( t \right) = 0$ or using approximation methods.

Answer 2

Assuming $x(t)=0\ \forall t\lt d$, and using the Laplace transform, we have

$$ \hat x(s) = \frac{x_0}{(s+h-r e^{-s d})} = \frac{x_0}{s+h}\frac{1}{1-\frac{r e^{-s d}}{(s+h)}} $$

now assuming $|r e^{-s d}| < |s+h|$ we have

$$ \hat x(s) = \frac{x_0}{s-h}\sum_{k=0}^{\infty}\left(\frac{r e^{-s d}}{s+h}\right)^k $$

and

$$ \mathcal{L}^{-1}\left[\left(\frac{r e^{-s d}}{s+h}\right)^k\right]=\frac{r^k (t-d k)^{k-1} e^{d h k-h t} \theta (t-d k)}{\Gamma (k)} $$

where $\theta(\cdot)$ is the Heaviside Theta function.

Follows a MATHEMATICA script computing an example for $r = 0.2, h = 0.3, x_0 = 1, d = 2$.

n = 5;
parms = {r -> 0.2, h -> 0.3, x0 -> 1, d -> 2};
xt = InverseLaplaceTransform[x0/(s + h) Sum[(r Exp[-s d]/(s + h))^k, {k, 0, n}], s, t];
xt0 = xt /. parms;
Plot[xt0, {t, 0, n d} /. parms, PlotStyle -> Blue, PlotRange -> All]

For $r=0.2, h=0.3$

For $r=0.2, h = 0.2$

For $r = 0.2, h = 0.1$