Analytic solutions of linear neutral and non-neutral delay differential equations using the Laplace transform method: featuring higher order poles and resonance

In this article, we extend the Laplace transform method to obtain analytic solutions for linear RDDEs and NDDEs which have real and complex poles of higher order. Furthermore, we present first-order linear DDEs that feature resonance phenomena. The procedure is similar to the one where all of the poles are order one, but requires one to use the appropriate modifications when using Cauchy’s residue theorem for the poles of higher order. The process for obtaining the solution relies on computing the relevant infinite sequence of poles and then determining the Laplace inverse, via the Cauchy residue theorem. For RDDEs, the poles can be obtained in terms of the Lambert W function, but for NDDEs,the complex poles, in most cases, must be computed numerically. We found that an important feature of first-order linear RDDES and NDDES with poles of higher order is that it is possible to incite the resonance phenomena, which in the counterpart ordinary differential equation cannot occur. We show that despite the presence of higher order poles or resonance phenomena, the solutions generated by the Laplace transform method for linear RDDEs and NDDEs that have higher order poles are still accurate.


Introduction
Delay differential equations (DDEs) appear in many applications from different fields of science [1][2][3][4][5][6][7][8][9][10]. For instance, several studies of climate models related to the global works have obtained resonance for nonlinear DDEs but not for linear ones [6,11,41,42]. We show that despite the presence of higher order real and complex poles or resonance phenomena, the solutions generated by the LTM for linear RDDEs and NDDEs that have higher order poles are accurate and especially for larger t values. Moreover, the steady-state behavior of the solutions can be calculated by the LTM. On the other hand, oftentimes the MoS can just produce the solution for a short time.
The remainder of this article proceeds as follows. In the next section, we develop the procedure to obtain the solutions of linear RDDEs and NDDEs with higher order poles by means of the LTM. In Sect. 3, we present numerical results for linear RDDEs and NDDEs with higher order poles. We measure the accuracy of the solutions produced by the LTM using the solutions obtained by the MoS. In Sect. 4, we present our conclusions.

Solving linear DDEs considering the presence of higher order poles
Here, we develop the procedure to obtain the solutions of linear RDDEs and NDDEs with higher order poles by means of the LTM. The process to obtain the solution by means of the LTM can be summarized in a series of steps; (i) first, the LT is applied to the DDE, (ii) we compute all the real poles and a finite number of complex poles of the denominator of the solution in the s-space, and (iii) we compute the ILT of the solution in the s-space by using Cauchy's residue theorem [39,40]. The special case of poles of order M = 2 is briefly discussed in [25], but only for the RDDE case. For the interested readers, we refer them to previous works [25,28,29], where the details of the methodology are presented in greater detail. However, we present some crucial steps in this section and in particular for the cases of poles of higher order for both RDDE and NDDE systems.
First, let us describe the role of the Lambert W function in solving linear DDEs. The characteristic equation of a linear DDE can take on the form of a quasi-polynomial which could have an infinite number of roots. The Lambert W Function provides a means for expressing these roots in terms of Lambert W branches. The equation takes on the form z = W (z)e W (z) where z is some complex number and W (z) is the complex multi-valued function solution of the equation. For each integer k, we denote W k (z) as the k-th branch of the Lambert function. The branch W 0 (z) is called the principal branch. For real numbers, the branches W 0 and W −1 can be used to solve the equation ye y = x with real numbers x and y if x ≥ − 1 e . The solution, y, can be found using the following conditions: We note that the second condition implies that in the interval − 1 e ≤ x < 0, the equation ye y = x has two real solutions. Further details and applications of the Lambert function can be found in literature such as [43,44]. Now, let us assume that f (z) is an analytic function in the complex region R defined by 0 < |z − z 0 | < d, or the neighborhood of a point z = z 0 ∈ C. An isolated singularity is a pole if f (z) can be written as [45,46]: where M ∈ N is the order of the pole, ψ(z 0 ) = 0 and ψ(z) is an analytic function in R.
Residues of Poles: From Eq. (1), if M = 1, then f (z) has a simple pole at z 0 . At a simple pole z 0 ∈ C, the residue of f (z) is given by Now, suppose f (z) can be written as f (z) = N (z) D(z) , where N (z) and D(z) are analytic functions in R and D(z) has a zero at z 0 . Let us assume that D(z) = (z − z 0 ) MD (z) whereD(z) is analytic in R andD(z 0 ) is nonzero. Using the Taylor series expansion of N (z) andD(z), the residue is given as [45,46]: Now if we have that M ≥ 2, it can be shown through the Taylor Series Expansion of ψ(z) about z 0 that the residue is given by [45,46]: After these basic definitions,we can proceed to solve linear RDDEs and NDDEs with the LTM.

Solving linear RDDEs in the presence of higher order poles
Let us consider the following linear RDDE with a non-homogeneous term g(t), Let us define L (y(t)) = Y (s). Considering the homogeneous case and employing the LT to Eq. (5), we obtain the LT transform, Y (s) = N (s) D(s) of y(t): Setting D(s) = 0, we find that the poles of Y (s) occur at s k = a + 1 τ W k τ be −aτ for k ∈ Z. From the definition of the Lambert function it can be shown that if τ be −aτ = − 1 e , then the poles s k are of order M = 1 and employing the ILT, one gets where the coefficients C −1,k 's are the calculated residues and the real part is taken due to the definition of the ILT. Now, if τ be −aτ = − 1 e , then we encounter a case where the real pole s = a − 1 τ is of order M = 2 (since s 0 = s −1 by definition of the W 0 and W −1 Lambert branches). Then, from Eq. (4), the residue for a pole s 0 = s −1 of order M = 2, we get Employing the ILT, one obtains Re lim for s 0 = a − 1 τ . For a linear non-homogeneous RDDE, one gets where G(s) is the LT of the function g(t). We denote the additional poles introduced by G(s) as s v,i for i = 1, 2, . . . , V , V ∈ Z + , with corresponding residues c v,i . Employing the ILT and the residue theorem, and assuming s k = s v,i for all k, i, we obtain the solution Re lim where It is possible that the roots of D(s) and the poles of G(s) may be the same. In which case, the term, Res{Y (s)e st ; s v,i }, is computed according to Eq. (4) for the poles that are the same (when s k = s v,i for some k, i). In this paper,we consider the case when a small subset of the poles of G(s) coincide with any of the infinite sequence of the roots of D(s).

Solving linear NDDEs in the presence of higher order poles
Let us consider the following linear NDDE, , of y(t): where . Therefore when c > 0, there are two real poles, s −1 = a and s 0 = ln(c)/τ. However, when c < 0, there will be only one real pole at s 0 = a. The second term in the factorization also provides one with the exact location of the complex poles above the x axis: and for k ∈ N. Depending on the LT of the non-homogeneous term, g(t), we may also obtain a repeated real or complex pole. However, we will see in the next section that a special case arises where we obtain a repeated real pole when ac + b = 0 and g(t) = 0. As in the RDDE case, we may encounter s k = s v,i for some k, i. In which case, the term, Res{Y (s)e st ; s v,i }, is computed according to Eq. (4). Otherwise, employing the ILT and the residue theorem, assuming s k = s v,i for all k, i, one gets where If r = s k or r = s v,i is a complex pole of order M = 2, then we have a complex conjugate pole and we could take advantage of this fact to compute the residue via, In this way, we reduce the number of operations and therefore the computation time, even though it is not significant. Notice that the procedure that takes longer is the development of the analytical solution which varies depending on the order of the poles, the non-homogeneous term, and the history function.

Results
In this section, we present the solutions of a variety of linear RDDEs and NDDEs with higher order poles. We assess the accuracy of the LTM solutions by comparing them with the analytical solutions provided by the MoS. We rely on Maple to obtain the solution by the LTM. For this section, let us introduce the next notation:

Linear RDDE examples
In this subsection, we consider several RDDEs in order to show different phenomena that can be obtained depending on the type of poles and the order of them.

Example 1 Consider the next RDDE:
where τ = 2. Taking the LT, we obtain the quasi-polynomial D(s) = (s − 1)e 2s + e 2 , which has a real pole at s 0 = 1 2 . We notice that the condition τ be −aτ = −e −1 is satisfied for the pole of order M = 2 given by s 0 = s −1 = a − 1 τ = 1 2 . The remaining poles and residues are as follows: Employing the ILT, one gets the LTM solution of y(t): The solution of Example 1 via MoS (over m = 5 intervals) and LTM (n = 100) is shown on the left-hand side of Fig. 1. We notice that as t → ∞, y(t) → −∞. The solution is unstable due to the presence of a positive real pole. The right-hand side of Fig. 1 shows the graph of the error. It can be observed that the solution generated by the LTM is accurate and its error decreases as t increases.
Example 2 Consider the next RDDE: where τ = 1. Taking the LT, we obtain the quasi-polynomial D(s) = (s− 1 3 )e s +e −2/3 , which has a real pole at s 0 = − 2 3 . We notice that the condition τ be −aτ = −e −1 is satisfied for the pole of order M = 2 given by s 0 = s −1 = a − 1 τ = − 2 3 . The remaining poles and residues are: Employing the ILT, one gets the LTM solution of y(t): The solution of Example 2 via MoS (over m = 10 intervals) and LTM (n = 100) is shown on the left-hand side of Fig. 2 . We notice that as t → ∞, y(t) → 0. The solution is stable since all of the poles have negative real parts. The right-hand side of Fig. 2 shows the graph of the error. It can be observed that the solution generated by the LTM is accurate and its error decreases as t increases.
where τ = 3. Denote G(s) = 1 2(s+1/6) as the LT of the function g(t) = e −t/6 2 . The additional pole introduced by G(s) is denoted as s v,1 = − 1 6 with residue c v,1 = 1 2 . Taking the LT, we obtain the quasi-polynomial D(s) = s − 1 3 e 3s + 1 3 , which has a real pole at s 0 = 0. We notice that the condition τ be −aτ = −e −1 is satisfied for the pole of order M = 2 given by s 0 = s −1 = a − 1 τ = 0. The remaining poles and residues are as follows: . Employing the ILT, one gets the LTM solution of y(t): The solution of Example 3 via MoS (over m = 5 intervals) and LTM (n = 100) is shown on the left-hand side of Fig. 3. We notice that as t → ∞, y(t) → −∞, and thus the system is unstable. The right-hand side of Fig. 3 shows the graph of the error. It can be observed that the solution generated by the LTM is accurate and its error decreases as t increases.
Example 4 Consider the next RDDE: where τ = 2. The characteristic equation . We also find that Re(s k ) < 0 for all integers k and thus, we expect the solution to be asymptotically stable. Employing the ILT, one gets the LTM solution of y(t): The solution of RDDE Example 4 via MoS (over m = 9 intervals) and LTM (n = 100 terms) in Maple is shown on the left-hand side of Fig. 4. We notice that as t → ∞, y(t) → 0. The solution is stable since all of the poles have negative real parts. The right-hand side of Fig. 4 shows the graph of the error between MoS and LTM. It can be observed that the solution generated by the LTM is accurate and its error decreases as t increases.
where τ = 3 4 . The characteristic equation of the system is: There are no real poles and there are an infinite number of complex poles, including s 0 = π i and s −1 = −π i. The LT of the non-homogeneous term, L{sin(π t)} = G(s) = π s 2 +π 2 , which means s 0 and s −1 are actually complex poles of order M = 2 with zero real part. The exact location of the remaining complex poles is: . We also find that Re(s k ) ≤ 0 for all integers k. However, due to the repeated complex pole with zero real part, we expect that the solution will oscillate with larger amplitudes as t increases. Employing the ILT, one gets the LTM solution of y(t):

The solution of RDDE Example 5 via MoS and LTM (n = 50 terms) in Maple is
shown on the left-hand side of Fig. 5 . We note that as t → ∞, where the coefficients A 1 , A 0 , B 1 and B 0 are constants. These terms are introduced because of the order 2 poles at s = ±iπ. Therefore, mathematically the solution assumes the same form as that of a massspring system, when the external force is chosen to induce resonance. The complex poles at s = ±iπ effectively constitute a natural frequency for the DDE (19). The resonance is caused by the inclusion of the non-homogeneous forcing function, with this same frequency. However, from the right-hand side of Fig. 5, it can be observed that the error of the LTM solution decreases as t → ∞. Note that despite the fact that the real part of all the poles is non-positive, the solution is unstable due to the repeated complex poles.
Since the condition ac + b = 0 is satisfied, then the exact location of the remaining complex poles is given by: s k = ln(c)+2kπ i τ = kπ i for k ∈ Z\{0}. Note, when k = −2, 2, then s −2 = s v,2 and s 2 = s v,1 are complex poles of order M = 2. Employing the ILT, one gets the LTM solution of y(t) (note, here N (s) = 0): The solution of NDDE Example 1 via MoS (over m = 8 intervals) and LTM (n = 50 terms) in Maple is shown on the left-hand side of Fig. 6. We notice that as t → ∞, y(t) → ∞. The right-hand side of Fig. 6 plots the error between MoS and LTM. Note that despite all poles have negative or zero real parts, the solution is unstable due to the resonance phenomenon. It can be observed that the solution generated by the LTM is accurate and its error decreases as t increases.
Since the condition ac + b = 0 is satisfied, then the exact location of the remaining complex poles is given by: s k = ln(c)+2kπ i τ = − ln(2) + 2kπ i for k ∈ Z\{0}. Note that Re(r 0 ) < 0, Re(r 1 ) < 0, and Re(s k ) < 0 for all k. Thus, we expect the solution to be asymptotically stable. Employing the ILT, one gets the LTM solution of y(t): The solution of NDDE Example 2 via MoS (over m = 20 intervals) and LTM (n = 100 terms) in Maple is shown on the left-hand side of Fig. 7. We notice that as t → ∞, y(t) → 0. This is due to the fact that all poles have negative real parts. The right-hand side of Fig. 7 plots the error between MoS and LTM.
Example 3 Consider the following NDDE: where τ = 1. The characteristic equation of the system is: which has a real pole at r 0 = 0 of order M = 2. In [28], it was shown that when c > 0 and for sufficiently large k ∈ Z\{0}, the remaining complex poles lie relatively close to the imaginary axis, and can be approximated with the initial guesses s k ≈ ln(c)+2kπ i τ ≈ −0.1054 + 2kπ i. Employing the ILT, one gets the LTM solution of y(t): The solution of NDDE Example 3 via MoS (over m = 24 intervals) and the LTM (n = 100 terms) in Maple is shown on the left of Fig. 8 . The terms in the infinite sum become negligible as t → ∞ because all of the complex poles have negative real parts. Therefore, as t → ∞, y(t) → ∞. More specifically Notice that the real pole at r 0 = 0, of order 2, introduces the linear growth term, resulting in a behavior that is in some way analogous to the resonance phenomenon. Figure 8b plots the error between MoS and LTM.

Conclusions
In this work, we extended the Laplace transform method to obtain the analytic solutions for linear RDDEs and NDDEs which have real and complex poles of higher order. The procedure is similar to the one where all of the poles are order one, but requires one to use the appropriate modifications when using Cauchy's Residue Theorem for the poles of higher order. This, as a consequence, introduces at least one additional term in the series solution which is typically of a different form to those at the simple poles. The process for obtaining the solution relies on computing the relevant infinite sequence of poles and then determining the Laplace inverse, via the Cauchy residue theorem. For RDDEs, the poles can be obtained in terms of the Lambert W function, but for NDDEs,the complex poles, in almost all cases, must be computed numerically. In [28],the authors derived a formula which provides an approximation for the location of the complex poles. In some cases, however, this formula does provide one with the exact location of the poles, and therefore, numerical methods are not required. We found that an important feature of first-order linear RDDES and NDDES with poles of higher order is that it is possible to incite the resonance phenomena, which, in the counterpart ordinary differential equation, cannot appear. This phenomenon can occur for instance when a complex pole has order two. For any NDDE which requires one to compute the complex poles numerically, it is essentially impossible to introduce a non-homogeneous term that guarantees a repeated complex pole. Since the complex poles are computed with numerical methods, almost always these complex poles have infinitely many digits, and therefore, it is difficult to have a non-homogeneous term that generates this type of complex pole. However, it is possible to select parameter sets that ensure that the formula for the approximate location of the complex poles provides one with the exact location of these poles. Therefore, we were able to design NDDEs with complex poles of higher order. We show that despite the presence of higher order real and complex poles or resonance phenomena, the LTM is able to generate accurate solutions of linear RDDEs and NDDEs. These solutions are even more accurate for larger t values. This is a main strength of the method since oftentimes the MoS can only provide the solution for the first few time intervals.
Author Contributions MS, GK,and GGP wrote the main manuscript text and prepared all figures. All authors reviewed the manuscript.
Funding Open access funding provided by SCELC, Statewide California Electronic Library Consortium

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.