Finding the analytical solution of a first order system with pure time delay

Question

I have a simple system and I am searching for the solution for f(t):

$$\frac{\partial f(t)}{\partial t} = c_1 \left( f(t) + g(t) + c_2 \right)$$.

It turns out that, in this system $g$ can be related to $f$ by a pure time delay:

$$g(t) = f(t-a)u(t-a)$$

where $u(t)$ is the Heaviside function.

Applying the Laplace transform to the equation above and isolating $F(s)$ leads to:

$$F(s) = \frac{f(0) + c_1c_2}{s - c_1(1 - e^{-as})}$$.

Since there is a complex exponential at the denominator I cannot use the inverse Laplace transform nor break this into partial fractions.

How can one find the solution of such an equation?

If I cannot directly resolve the equation, can I infer a given solution (from the numerical integration of this system) and identify the coefficients?

If the analytical solution does not exists, is there a proof that explains why?

Answer 1

Applying the Laplace transform to

$$\frac{\partial f(t)}{\partial t} = c_1 \left( f(t) + f(t-a)u(t-a) + c_2 \right)$$.

and assuming $f(t)=0,\ t < a$, we got

$$ F(s) = -\frac{((c_1c_2+s f_0)e^{as})}{s(c_1+(c_1-s)e^{as})}=-\frac{c_1c_2+s f_0}{s(c_1e^{-a s}+c_1-s)}=\frac{c_1c_2+s f_0}{s(s-c_1)(1-\frac{c_1e^{-as}}{s-c_1})} $$

and now making

$$ \frac{1}{(1-\frac{c_1e^{-as}}{c_1-s})}= \sum_{k=0}^{\infty} \left(\frac{c_1e^{-as}}{s-c_1}\right)^k $$

we can invert $F(s)\to f(t)$.

Note that the series is convergent for $|c_1e^{-as}|<|s-c_1|$ and

$$ \frac{c_1c_2+s f_0}{s(s-c_1)}=\frac{c_2+f_0}{s-c_1}-\frac{c_2}{s} $$

and

$$ \mathcal{L}^{-1}\left[\frac{\left(\frac{c_1e^{-as}}{s-c_1}\right)^k}{s-c_1}\right]=\frac{c_1^k \left(\left(\left(1-\frac{t}{a k}\right)^k-1\right)(-a k)^k\right) e^{c_1 (t-a k)} u (t-a k)}{k!} $$

Follows a MATHEMATICA script showing the approximation accuracy for n = 4 in blue and the precise integration in dashed red

n = 4;
parms = {c1 -> -1, c2 -> 1, a -> 2, f0 -> 0};
sum = Sum[((c1 Exp[-a s])/(s - c1))^k, {k, 0, n}];
tmax = (n+1) a /. parms;
Fs = (c1 c2 + s f0)/(s (s - c1)) sum /. parms;
ft = InverseLaplaceTransform[Fs, s, t];
gr1 = Plot[ft, {t, 0, tmax}, PlotRange -> All];
ode = f'[t] - c1 (f[t] + f[t - a] + c2);
ode0 = {ode == 0, f[0] == f0} /. parms;
ft = Quiet@NDSolve[ode0, f, {t, 0, tmax}][[1]];
gr2 = Plot[Evaluate[f[t] /. ft], {t, 0, tmax}, PlotStyle -> {Dashed, Red}, PlotRange -> All];
Show[gr1, gr2]

Answer 2

Solving

The differential equation is a linear DDE of order one with a non-constant coefficient. We can get rid of this one by splitting the DDE into two intervals: \begin{align*} \frac{\operatorname{d}f\left( t \right)}{\operatorname{d}t} &= c_{1} \cdot \left( f\left( t \right) + f\left( t - a \right) \cdot \theta\left( t - a \right) + c_{2} \right)\\ \frac{\operatorname{d}f\left( t \right)}{\operatorname{d}t} &= c_{1} \cdot f\left( t \right) + c_{1} \cdot f\left( t - a \right) \cdot \theta\left( t - a \right) + c_{1} \cdot c_{2}\\ \frac{\operatorname{d}f\left( t \right)}{\operatorname{d}t} &= \begin{cases} c_{1} \cdot f\left( t \right) + c_{1} \cdot c_{2}, &\text{for } t < a\\ c_{1} \cdot f\left( t \right) + c_{1} \cdot f\left( t - a \right) + c_{1} \cdot c_{2}, &\text{for } t \geq a\\ \end{cases}\\ \end{align*}

Now there are two linear DEs of order $1$ but constant coefficients. We can trivially for $t < a$ (because it's an order $1$ ODE) where $c$ can be any constant: \begin{align*} f\left( t \right) &= c \cdot e^{c_{1} \cdot t} - c_{2}\\ \end{align*}

For $t \geq a$ we could say $f\left( t \right) = g\left( t \right) + b$. We would get: \begin{align*} \frac{\operatorname{d}g\left( t \right)}{\operatorname{d}t} &= c_{1} \cdot g\left( t \right) + c_{1} \cdot g\left( t - a \right) + 2 \cdot c_{1} \cdot b + c_{1} \cdot c_{2}\\ \end{align*}

If $b = -\tfrac{1}{2} \cdot c_{2}$ we get a homogeneous DDE: \begin{align*} \frac{\operatorname{d}g\left( t \right)}{\operatorname{d}t} &= c_{1} \cdot g\left( t \right) + c_{1} \cdot g\left( t - a \right)\\ \end{align*}

Substituting $g\left( t \right) = e^{\lambda \cdot t}$ we get: \begin{align*} \frac{\operatorname{d}e^{\lambda \cdot t}}{\operatorname{d}t} &= c_{1} \cdot e^{\lambda \cdot t} + c_{1} \cdot e^{\lambda \cdot \left( t - a \right)}\\ \lambda \cdot e^{\lambda \cdot t} &= c_{1} \cdot e^{\lambda \cdot t} + c_{1} \cdot e^{\lambda \cdot t} \cdot e^{-\lambda \cdot a}\\ \lambda &= c_{1} + c_{1} \cdot e^{-\lambda \cdot a}\\ \end{align*}

This is an classical equation wich can be solved via the Lambert W-function: \begin{align*} \lambda &= c_{1} + c_{1} \cdot e^{-a \cdot \lambda} \tag{$-c_{1}$}\\ \lambda - c_{1} &= c_{1} \cdot e^{-a \cdot \lambda}\\ \lambda - c_{1} &= c_{1} \cdot e^{-a \cdot c_{1} - a \cdot \left( \lambda - c_{1} \right)} \tag{$\cdot e^{a \cdot \left( \lambda - c_{1} \right)}$}\\ \left( \lambda - c_{1} \right) \cdot e^{a \cdot \left( \lambda - c_{1} \right)} &= c_{1} \cdot e^{-a \cdot c_{1}} \tag{$\cdot a$}\\ a \cdot \left( \lambda - c_{1} \right) \cdot e^{a \cdot \left( \lambda - c_{1} \right)} &= a \cdot c_{1} \cdot e^{-a \cdot c_{1}} \tag{$\operatorname{W}_{k}\left( \cdot \right)$}\\ a \cdot \left( \lambda - c_{1} \right) &= \operatorname{W}_{k}\left( a \cdot c_{1} \cdot e^{-a \cdot c_{1}} \right) \tag{$\div a$ then $+ c_{1}$}\\ \lambda &= \frac{1}{a} \cdot \operatorname{W}_{k}\left( a \cdot c_{1} \cdot e^{-a \cdot c_{1}} \right) + c_{1}\\ \end{align*} where $k \in \mathbb{Z}$. So:

\begin{align*} f\left( t \right) = \begin{cases} c \cdot e^{c_{1} \cdot t} - c_{2}, &\text{for } t < a\\ e^{\left( \frac{1}{a} \cdot \operatorname{W}_{k}\left( a \cdot c_{1} \cdot e^{-a \cdot c_{1}} \right) + c_{1} \right) \cdot t} - \frac{c_{2}}{2}, &\text{for } t \geq a\\ \end{cases} \end{align*}

or $$\fbox{$\begin{align*} f\left( t \right) = \begin{cases} c \cdot e^{c_{1} \cdot t} - c_{2}, &\text{for } t < a\\ \sum\limits_{k = -\infty}^{\infty}\left[ d_{k} \cdot e^{\left( \frac{1}{a} \cdot \operatorname{W}_{k}\left( a \cdot c_{1} \cdot e^{-a \cdot c_{1}} \right) + c_{1} \right) \cdot t} \right] - \frac{c_{2}}{2}, &\text{for } t \geq a\\ \end{cases} \end{align*}$}$$ where $d_{k} \in \mathbb{C}$ and $k \in \mathbb{Z}$ (after all, the terms with the Lambert W-function are homogeneous solutions and there are infinitely many in both directions).

In Mathematica you could check the solution $g\left( t \right)$ and it's values like this:

g[t_] := Subscript[d, 1] Exp[(1/a LambertW[k, a Subscript[c, 1] Exp[-a Subscript[c, 1]]] + Subscript[c, 1]) t];
Subscript[c, 1]:= 1; Subscript[c, 2] := 20; Subscript[d, 1] := 2; k := 1; a := 1; x := 20;
N[g'[x] == Subscript[c, 1] g[x] + Subscript[c, 1] g[x - a]]

It gives you the output True. For some reason this is not the case when adding the constant term in the DE an the functional equation even if it has to be so.