Comparison of Symbolic Computations for Solving Linear Delay Differential Equations Using the Laplace Transform Method

Abstract

In this paper, we focus on investigating the performance of the mathematical software program Maple and the programming language MATLAB when using these respective platforms to compute the method of steps (MoS) and the Laplace transform (LT) solutions for neutral and retarded linear delay differential equations (DDEs). We computed the analytical solutions that are obtained by using the Laplace transform method and the method of steps. The accuracy of the Laplace method solutions was determined (or assessed) by comparing them with those obtained by the method of steps. The Laplace transform method requires, among other mathematical tools, the use of the Cauchy residue theorem and the computation of an infinite series. Symbolic computation facilitates the whole process, providing solutions that would be unmanageable by hand. The results obtained here emphasize the fact that symbolic computation is a powerful tool for computing analytical solutions for linear delay differential equations. From a computational viewpoint, we found that the computation time is dependent on the complexity of the history function, the number of terms used in the LT solution, the number of intervals used in the MoS solution, and the parameters of the DDE. Finally, we found that, for linear non-neutral DDEs, MATLAB symbolic computations were faster than Maple. However, for linear neutral DDEs, which are often more complex to solve, Maple was faster. Regarding the accuracy of the LT solutions, Maple was, in a few cases, slightly better than MATLAB, but both were highly reliable.

1. Introduction

Often, engineers and scientists may need to choose a software that is specific to their application or modeling purpose. Some software languages may have different features, speed/computational complexity, tolerances, architecture, libraries, ease of use, etc., that can affect the accuracy of the model. Here, we investigated the performance of the mathematical software program Maple and the programming language MATLAB when applying the method of steps (MoS) and the Laplace transform (LT) methods to solve linear delay differential equations (DDEs). To make the comparison fair, we used the benefits offered by these two programs for symbolic computation, which allows us to derive the solutions of these DDEs in the ways we would think about carrying them out by hand with variables, symbols, functions, and other mathematical formulas. An advantage of using symbolic computation is that it eliminates cases of roundoff errors in intermediate calculations since parameter variables are carried throughout the calculations using infinite precision arithmetic [1,2,3]. However, in some cases, using symbolic computations can lead to slower execution times depending on the program used and the process [4,5,6]. Symbolic computation has been used and proposed for solving many different problems [7,8,9,10,11].

In this paper, we focus on comparing the symbolic mathematical computations of Maple and MATLAB to obtain analytical solutions of linear neutral and non-neutral DDEs by using the LT method. The accuracy of the LT solutions may be determined by comparing them with the MoS solutions [12,13,14,15,16,17]. However, the MoS has a very important limitation since it is a stepwise process that requires solving the relevant differential equation in each respective time interval, which can be computationally expensive. Moreover, for many types of DDEs, we will often encounter the swelling phenomena, i.e., the complexity of the closed-form solution increases rapidly in each successive time interval, making the computations unfeasible [14,16]. On the other hand, the LT method produces a solution that can be evaluated at any particular value of time. Thus, we can compute the approximated solution by evaluating it at any point. Other advantages of the analytical solutions, such as the long-term behavior of the solution, can be seen in [14,15,18,19,20,21,22]. The errors of the LT solutions usually decrease as time increases. On the other hand, for classical numerical solutions, the computation increases as $t\to \infty$ [14,23,24,25]. However, numerical schemes that solve linear and nonlinear DDEs are useful and, in some cases, are easy to implement. For instance, in [26], the authors developed a numerical scheme based on a fitted difference scheme with the use of exponential basis functions and interpolating quadrature rules to solve first order linear DDEs. Some researchers have proposed different collocation methods to solve linear and nonlinear DDEs [27,28]. For instance, in [29], the Haar wavelet collocation method was applied to obtain the numerical solution of a particular class of DDEs. The method is also applied to systems involving DDEs. In [30], the authors proposed a numerical method for solving linear fractional differential equations using discretization at the Jacobi–Gauss collocation points. Other methods that have been applied that involve more general DDEs can be found in [31,32]. In these works, the authors deal with the differential-difference equation in an infinite-dimensional space and with unbounded domains. Moreover, they consider mixed-type equations that include both forward and backward shifts of the argument. These forms are generated by different temporal and spatial discretizations. It is important to mention that one might expect that some of the techniques and results of retarded differential equations could be applied to mixed-type equations, but, given that the independent variable in this paper is strictly temporal, this may require one to introduce some modifications [32,33]. In addition, these last works deal with idealized nonlinearities in order to obtain a more tractable problem.

It is important to point out that, even for linear DDEs, it might be challenging to obtain an analytical solution. Moreover, in most cases, the neutral delay differential equations (NDDEs) are more difficult to solve due to the appearance of a time delay on the derivative of the state variable. Regarding solving ordinary differential equations (ODEs) and DDEs, the latter is usually more difficult to solve since the derivatives of the state variables depend on the values of the state variable from the past. The classic approach used to obtain the solutions of DDEs is the MoS, but, in many cases, it is unfeasible to advance further in time and obtain the solution in the long term [14,15,34,35,36]. In [14], the authors studied the accuracy of the LT for linear DDEs and, in particular, for NDDEs. In [20], the authors developed a method for solving linear and nonlinear DDEs with a proportional delay using a combination of the LT and the differential transformation method (DTM). The authors mentioned that the convergence rate can be improved by the combination of the two methods. The LT method has also been used to solve linear second order DDEs [37]. However, complex analytical computations and the use of the residue theorem are required [14,37]. Linear DDEs have different applications and have been applied in many problems [17,38,39].

In order to compare the symbolic computations involved in obtaining the solutions of neutral and non-neutral linear DDEs using the LT method, we selected Maple and MATLAB software since they are well-known in the scientific community [3,40]. In Maple, we made use of the symbolic functions diff, dsolve, fsolve, eval, and subs, to name a few, whereas, in MATLAB, we utilized functions from the symbolic math and optimization toolboxes such as fsolve, subs, vpa, and vpasolve. We explored several examples of DDEs of neutral and non-neutral type in order to compare the symbolic computation of Maple and MATLAB. For each example, we computed the solution over some finite time interval using the MoS briefly described in Section 2.1 and the LT method described in Section 3. First, we explored two examples of non-neutral DDEs and then we followed with several examples of neutral DDEs. For each example, we compared: (1) the solutions of the two methods from each program software, (2) the errors associated with each method and software, and (3) the computation times. DDEs have been used in many different fields, including astronomy and epidemiology [17,41,42,43]. For some applications, such as control or communication systems, the computation time may be a critical factor [44,45,46,47]. All computations were performed on a computer system with an Intel(R) Core(TM) i9-10900K CPU @ 3.70 GHz processor and 64 GB RAM.

The organization of this paper is as follows. In Section 2, preliminaries for DDEs, NDDEs, and the MoS for solving linear DDEs are presented. In Section 3, we briefly present the main aspects of the LT method with the relevant equations and formulas. In Section 4, we compare the symbolic computations to obtain the solution of linear DDEs using the LT method. We present several examples related to DDEs and their corresponding numerical solutions. In Section 5, conclusions regarding the symbolic computations required to obtain the solution of DDEs using the LT method are provided.

2. Preliminaries

In this section, we present some preliminary definitions that are needed to solve the different linear DDEs presented in this work. Some preliminaries include background concepts, the Lambert function, the residue theorem, and the MoS. Familiarity with the basic concepts of complex analysis is assumed.

2.1. Method of Steps (MoS)

The traditional method for solving DDEs is with an elementary method called the method of steps [17,34,35,48,49,50]. It is well-known for being a tedious process when being worked out by hand, but the complexity can be drastically simplified with symbolic toolboxes in computer programs [15]. Let us consider a NDDE with fixed delay $\tau$:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =f(t,y\left(t\right),y(t-\tau ),y^{\prime }(t-\tau )),\hfill & & t>t_0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right),\hfill & & t\in [-\tau +t_0,t_0],\hfill \\ \hfill y\left(t\right)& =y_0,\hfill & & t=t_0,\hfill \end{array}$$

(1)

where $t_0$ is the initial time value, $\tau$ is a constant/discrete delay, and $y_0$ is some constant. For DDE initial value problems (IVPs), we require a history function, ${\varphi }_0\left(t\right):[-\tau ,t_0]\to \mathbb{R}$, and an initial value, $y\left(t_0\right)=y_0$ for $t=t_0$, to solve the problem. This condition implies that DDEs are in fact infinite-dimensional problems since we define an infinite set of initial conditions between $t\in [-\tau ,t_0]$.

Using integration properties on the interval $[t_0,t_0+\tau ]$, the solution is then given by $y_1\left(t\right)$, which is a solution to the following NDDE:

$$\begin{array}{cccc}\hfill y_1^{\prime }\left(t\right)& =f(t,y_1\left(t\right),{\varphi }_0(t-\tau ),{\varphi }_0^{\prime }(t-\tau )),\hfill & & t\in [t_0,t_0+\tau ],\hfill \\ \hfill y_1\left(t\right)& ={\varphi }_0\left(t\right),\hfill & & t=t_0.\hfill \end{array}$$

Repeating this process for successive intervals and using the solution over the past interval, we find the following general derivative formula for the nth partition of the IVP ($n\in \mathbb{N}$):

$$\begin{array}{cccc}\hfill y_n^{\prime }\left(t\right)& =f(t,y_n\left(t\right),y_{n-1}(t-\tau ),y_{n-1}^{\prime }(t-\tau )),\hfill & & t\in \left[\right(n-1)\tau ,n\tau ],\hfill \\ \hfill y_n\left(t\right)& =y_{n-1}\left(t\right),\hfill & & t=(n-1)\tau .\hfill \end{array}$$

The solutions of ${\left\{y_n^{\prime }\left(t\right)\right\}}_{n=1}^{\infty }$ form a piecewise solution of the original NDDE (1) on the interval $t\in [0,n\tau )$.

Even though the MoS provides an analytical solution for DDEs, it has some drawbacks, including time consumption and the complexity of the integrals to be solved over each interval. This method also requires full knowledge of the history before the approach can be used. Since it is a successive process, it can be quite difficult to observe the behavior of the solution or to perform stability analysis of $y\left(t\right)$ as $t\to \infty$.

2.2. Lambert W Function

From ODEs, the roots of the characteristic equation provide information about the solution of a linear ODE or system of linear ODEs. Often times, this equation is a simple polynomial that is easy to solve with algebraic operations or use of the quadratic formula. A linear DDE or system of linear DDEs can be characterized in a similar manner to ODEs by looking at the characteristic equation. Take, for example, the following linear DDE with discrete time delay $\tau$:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =cy(t-\tau ),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right),\hfill & & t\in [-\tau ,0],\hfill \\ \hfill y\left(t\right)& =y_0,\hfill & & t=0,\hfill \end{array}$$

where $y_0$ is the initial value. The characteristic polynomial of the above DDE takes the form of a quasi-polynomial due to the introduction of the exponential, i.e., $P\left(\lambda \right)=\lambda e^{\tau \lambda }-c{y}_0$. By setting $P\left(\lambda \right)=0$, there could be an infinite number of poles. Through simplifying, we obtain $\lambda e^{\tau \lambda }=c{y}_0$, which takes the form of a Lambert W function with the equation form $z=W\left(z\right)e^{W\left(z\right)}$, where z is some complex number and $W\left(z\right)$ is the multi-valued function solution of the equation.

Lambert functions are used when the unknown of an equation appears as both a base and exponent. Equations such as $\lambda e^{\tau \lambda }=c{y}_0$ cannot be solved with algebraic, logarithmic, or exponential properties. For certain values of z, this equation can have an infinite number of solutions called branches of W, and are denoted by $W_k\left(z\right)$ for $k\in \mathbb{Z}$. One of the branches, called the principal branch, is analytic at zero and is denoted by $W_0$ ($k=0$). For real numbers, the branches $W_0$ and $W_{-1}$ can be used to solve the equation $y{e}^y=x$ with real numbers x and y if $x\ge -\frac{1}{e}$. The solution, y, can be found using the following conditions.

$\left\{\begin{array}{cc}y=W_0\left(x\right)\hfill & \mathrm{if}x\ge 0\hfill \\ y=W_{-1}\left(x\right)\mathrm{or}y=W_0\left(x\right)\hfill & \mathrm{if}-\frac{1}{e}\le x<0.\hfill \end{array}\right.$

We note that the second condition implies that, in the interval $-\frac{1}{e}\le x<0$, the equation $y{e}^y=x$ has two real solutions.

In terms of our given problem, the poles of the quasi-polynomial are given by ${\lambda }_k=\frac{1}{\tau }{W}_k\left(c{y}_0\tau \right)$. The analytical solution to the DDE, then, can be expressed in terms of an infinite number of branches, given as

$$y\left(t\right)=\sum _{k=-\infty }^{\infty }{e}^{{\lambda }_{k}t}{C}_k,{\lambda }_k=\frac{1}{\tau }{W}_k\left(c{y}_0\tau \right)$$

where the $C_k$s are the coefficients determined from the history function, ${\varphi }_0\left(t\right)$, and the initial value, $y_0$. One can analyze how changes in the parameters c, $y_0$, and $\tau$ affect the real part of the eigenvalues of the system and, in turn, how that impacts the stability of the solution. Further details and applications of the Lambert function can be found in literature such as [51,52].

2.3. Cauchy’s Residue Theorem

The LT method for solving DDEs requires the use of Cauchy’s residue theorem, which is a powerful tool discussed in complex analysis for evaluating complex real integrals or infinite series [53,54]. Here, it was used to find the solution of a DDE when converting from the s-space back to the t-domain.

Suppose $f\left(z\right)$ is an analytic function in the complex region $\mathcal{R}$ defined by $0<|z-z_0|<d$, or the neighborhood of a point $z=z_0\in \mathbb{C}$. Then, $z_0$ is an isolated singular point of $f\left(z\right)$. The function $f\left(z\right)$ can be represented by the Laurent series expansion as [54]:

$$\begin{array}{c}\hfill f\left(z\right)=\sum _{n=-\infty }^{\infty }{C}_n{(z-z_0)}^n\end{array}$$

where

$$\begin{array}{c}\hfill C_n=\frac{1}{2\pi i}{\oint }_C\frac{f\left(z\right)}{{(z-z_0)}^{n+1}}dz\end{array}$$

are the coefficients and C is a simple closed contour in $\mathcal{R}$. The principal part of the series is given by the negative portion, ${\displaystyle \sum _{n=-\infty }^{-1}}{C}_n{(z-z_0)}^n$, where $C_{-1}$ is called the residue of $f\left(z\right)$ at $z_0$. An isolated singularity is a pole if $f\left(z\right)$ is of the form [53,54]:

$$\begin{array}{c}\hfill f\left(z\right)=\frac{\psi \left(z\right)}{{(z-z_0)}^M}\end{array}$$

(2)

where $M\in \mathbb{N}$ is the order of the pole, $\psi \left(z\right)$ is an analytic function in $\mathcal{R}$, and $\psi \left(z_0\right)\ne 0$.

Residues of Poles: If $M=1$, then $f\left(z\right)$ has a simple pole at $z_0$. At a simple pole $z_0\in \mathbb{C}$, the residue of $f\left(z\right)$ is given by [54]:

$$\begin{array}{c}\hfill C_{-1}=\mathrm{Res}(f\left(z\right),z_0)=\underset{z\to z_0}{lim}\left((z-z_0)f\left(z\right)\right)=\psi \left(z_0\right).\end{array}$$

(3)

If the limit does not exist, then there exists an essential singularity at $z_0$. If the limit is zero, then $f\left(z\right)$ is analytic at $z_0$, or $z_0$ is a removable singularity. If the limit is infinity, then the order is greater than $M=1$ [53,54].

Now, suppose $f\left(z\right)$ can be expressed as a rational function, $f\left(z\right)=\frac{N\left(z\right)}{D\left(z\right)}$, where $N\left(z\right)$ and $D\left(z\right)$ are analytic functions in $\mathcal{R}$ and $D\left(z\right)$ has a zero at $z_0$. Let $D\left(z\right)={(z-z_0)}^M\widehat{D}\left(z\right)$, where $\widehat{D}\left(z_0\right)$ is analytic in $\mathcal{R}$ and $\widehat{D}\left(z_0\right)$ is nonzero. Using the Taylor series expansion of $N\left(z\right)$ and $\widehat{D}\left(z\right)$, the residue is given as [53,54]:

$$\begin{array}{cc}\hfill C_{-1}& =\psi \left(z_0\right)=\frac{N\left(z_0\right)}{D^{\prime }\left(z_0\right)}.\hfill \end{array}$$

(4)

If $M\ge 2$, then $f\left(z\right)$ has a pole of order M at $z_0$. From Equation (2), the residue of $f\left(z\right)$ at $z_0$ is given by [53,54]:

$$\begin{array}{cc}\hfill C_{-1}& =\frac{1}{(M-1)!}\underset{z\to z_0}{lim}\left(\frac{d^{M-1}}{d{z}^{M-1}}\left({(z-z_0)}^{M}f\left(z\right)\right)\right)=\frac{{\psi }^{(M-1)}\left(z_0\right)}{(M-1)!}.\hfill \end{array}$$

Residue Theorem: Suppose C is a simple closed contour and let $f\left(z\right)$ be an analytic function inside and on C, excluding a finite number (K) of isolated singularities inside C. Then,

$$\begin{array}{c}\hfill {\oint }_{C}f\left(z\right)dz=2\pi i\sum _{n=0}^{K-1}\mathrm{Res}(f\left(z\right),z_n).\end{array}$$

(5)

The use of Cauchy’s residue theorem in the LT method comes from the definition of the inverse Laplace transform (ILT):

$$\begin{array}{c}\hfill f\left(t\right)={\mathcal{L}}^{-1}\left\{F\left(s\right)\right\}=\frac{1}{2\pi i}\underset{c-i\infty }{\overset{c+i\infty }{\int }}F\left(s\right)e^{st}ds\end{array}$$

(6)

where the integration is over a vertical path in the complex plane defined at $c=$ Re$\left(s\right)$. From complex analysis, the function $e^{st}$ is entire and, thus, we require the vertical path to be to the right of all of the poles of $F\left(s\right)$. In other words, the real part of all of the poles needs to be to the left of c. Let us denote this vertical path $\gamma$ shown in Figure 1. To use the residue theorem to evaluate the above integral, we construct a contour path, $C_R$, enclosing all of the poles with radius R as shown in Figure 1 going in a CCW direction. From this contour and Cauchy’s residue theorem, we obtain

$$\begin{array}{cc}\hfill {\oint }_{C}F\left(s\right)e^{st}ds& ={\int }_{\gamma }F\left(s\right)e^{st}ds+{\int }_{C_R}F\left(s\right)e^{st}ds\hfill \\ \hfill & =\underset{c-iR}{\overset{c+iR}{\int }}F\left(s\right)e^{st}ds+\underset{C_R}{\int }F\left(s\right)e^{st}ds=2\pi i\sum _{n=0}^{N-1}\mathrm{Res}(F\left(s\right)e^{st};s_n),\hfill \end{array}$$

which implies that

$$\begin{array}{cc}\hfill f\left(t\right)& =\frac{1}{2\pi i}\underset{R\to \infty ,N\to \infty }{lim}\left(\underset{c-iR}{\overset{c+iR}{\int }}F\left(s\right)e^{st}ds\right)=\sum _{n=0}^{\infty }\mathrm{Res}(F\left(s\right)e^{st};s_n)-\frac{1}{2\pi i}\underset{R\to \infty }{lim}\underset{C_R}{\int }F\left(s\right)e^{st}ds.\hfill \end{array}$$

(7)

Let $I_{C_R}=\underset{C_R}{\int }F\left(s\right)e^{st}ds$. We need to show that $\underset{R\to \infty }{lim}\left(I_{C_R}\right)=\underset{R\to \infty }{lim}\underset{C_R}{\int }F\left(s\right)e^{st}ds=0$. Achieving this will give:

$$\begin{array}{c}\hfill f\left(t\right)=\frac{1}{2\pi i}\underset{c-i\infty }{\overset{c+i\infty }{\int }}F\left(s\right)e^{st}ds=\sum _{n=0}^{\infty }\mathrm{Res}(F\left(s\right)e^{st};s_n)\end{array}$$

where c is a fixed real constant s.t. $F\left(s\right)$ has a finite number of poles located to the left of the line $\mathrm{Re}\left(s\right)=c$.

Let $F\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}$. For $t>0$, we enclosed in the left half plane where $s=c+R{e}^{i\theta }$,$\theta \in [\pi /2,3\pi /2]$. Then, $ds=iR{e}^{i\theta }d\theta$, and we obtain

$$\begin{array}{cc}\hfill I_{C_R}& =\underset{\pi /2}{\overset{3\pi /2}{\int }}{e}^{(c+R{e}^{i\theta })t}·F(c+R{e}^{i\theta })·\left(iR{e}^{i\theta }\right)d\theta \hfill \\ \hfill \Rightarrow \left|I_{C_R}\right|& \le \underset{\pi /2}{\overset{3\pi /2}{\int }}\left|e^{ct+Rt{e}^{i\theta }}\right|·\left|F(c+R{e}^{i\theta })\right|·\left|iR{e}^{i\theta }\right|d\theta \hfill \\ \hfill & =\underset{\pi /2}{\overset{3\pi /2}{\int }}\left|e^{ct}{e}^{Rt\left(\cos \right(\theta )+i\sin (\theta \left)\right)}\right|·\left|F(c+R{e}^{i\theta })\right|·\left|i(s-c)\right|d\theta \hfill \end{array}$$

We find that

$$\begin{array}{c}\hfill \left|e^{ct}{e}^{Rt\left(\cos \right(\theta )+i\sin (\theta \left)\right)}\right|=e^{ct}{e}^{Rt\cos \left(\theta \right)}\end{array}$$

and

$$\begin{array}{c}\hfill \left|i\right(s-c\left)\right|=\left|i\right|·|s-c|\le \left|s\right|+\left|c\right|<2\left|s\right|\end{array}$$

for sufficiently large $\left|s\right|$. Then, we have

$$\begin{array}{cc}\hfill \left|I_{C_R}\right|& \le 2e^{ct}\underset{\pi /2}{\overset{3\pi /2}{\int }}{e}^{Rt\cos \left(\theta \right)}·\left|sF\left(s\right)\right|d\theta \hfill \end{array}$$

In the integrand, we note that $\cos \left(\theta \right)$ is negative for $\theta \in (\pi /2,3\pi /2)$. Then, as $R\to \infty$, we can observe that $|I_{C_R}|\to 0$, which implies that $\underset{R\to \infty }{lim}\left(I_{C_R}\right)=\underset{R\to \infty }{lim}\underset{C_R}{\int }F\left(s\right)e^{st}ds=0$. Then, we have shown that the ILT of $F\left(s\right)$ can be represented using residues and is given by

$$\begin{array}{cc}\hfill f\left(t\right)& =\sum _{n=0}^{\infty }\mathrm{Res}(F\left(s\right)e^{st};s_n).\hfill \end{array}$$

(8)

3. Laplace Transform (LT) Method

We denote the LT of the solution $y\left(t\right)$ as $\mathcal{L}\left\{y\right(t\left)\right\}=Y\left(s\right)$. Then, for a given delay term in the state variable, the LT is given by (such that $t_0=0$)

$$\begin{array}{cc}\hfill \mathcal{L}\left\{y\right(t-\tau \left)\right\}& =\underset{0}{\overset{\infty }{\int }}y(t-\tau )e^{-st}dt=e^{-s\tau }\underset{-\tau }{\overset{\infty }{\int }}y\left(v\right)e^{-sv}dv\hfill \\ \hfill & =e^{-s\tau }Y\left(s\right)+e^{-s\tau }\underset{-\tau }{\overset{0}{\int }}y\left(v\right)e^{-sv}dv.\hfill \end{array}$$

(9)

When the delay is present in the derivative of the state variable, the LT is given by (such that $t_0=0$):

$$\begin{array}{cc}\hfill \mathcal{L}\left\{y^{\prime }(t-\tau )\right\}& =\underset{0}{\overset{\infty }{\int }}{y}^{\prime }(t-\tau )e^{-st}dt=e^{-s\tau }\underset{-\tau }{\overset{\infty }{\int }}{y}^{\prime }\left(v\right)e^{-sv}dv\hfill \\ \hfill & =s{e}^{-s\tau }Y\left(s\right)+e^{-s\tau }\left[\underset{-\tau }{\overset{0}{\int }}{y}^{\prime }\left(v\right)e^{-sv}dv-y\left(0\right)\right]\hfill \\ \hfill & =s{e}^{-s\tau }Y\left(s\right)-y(-\tau )+s{e}^{-s\tau }\underset{-\tau }{\overset{0}{\int }}y\left(v\right)e^{-sv}dv.\hfill \end{array}$$

(10)

The LT method uses Equations (9) and (10) as a part of the process to obtain the solutions of the linear NDDEs. For further details, interested readers are referred to [14,15].

3.1. Theory for Linear Non-Neutral DDE Case

First, we consider a simple linear non-neutral DDE case as shown below, where the solution is derived using the ILT and the residue formula as described in Section 2.3:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =ay\left(t\right)+by(t-\tau )+g\left(t\right),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right),\hfill & & t\in [-\tau ,0].\hfill \end{array}$$

(11)

Considering the homogeneous case, and applying the LT to Equation (11) (using Equation (9)), we obtain the LT transform, $Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}$, of $y\left(t\right)$:

$$Y\left(s\right)=\frac{y\left(0\right)e^{s\tau }+b\underset{-\tau }{\overset{0}{\int }}{\varphi }_0\left(v\right)e^{-sv}dv}{(s-a)e^{s\tau }-b}.$$

(12)

Setting $D\left(s\right)=0$, we find that the poles of $Y\left(s\right)$ occur at $s_k=a+\frac{1}{\tau }{W}_k\left(\tau b{e}^{-a\tau }\right)$ for $k\in \mathbb{Z}$. Recall from the definition of the Lambert function that, if the argument $\tau b{e}^{-a\tau }=-\frac{1}{e}$, then we encounter a case where the real pole $s=a-\frac{1}{\tau }$ is of order $M=2$ (since $s_0=s_{-1}$ by definition of the $W_0$ and $W_{-1}$ Lambert branches). If $\tau b{e}^{-a\tau }\ne -\frac{1}{e}$, then the poles $s_k$ are of order $M=1$. By applying the ILT, for $M=1$, we obtain the solution

$$y\left(t\right)=\sum _{k\in \mathbb{Z}}\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right),$$

(13)

where the $c_k$s are the computed residues. Further details can be found in [14].

If the linear non-neutral DDE is non-homogeneous, we obtain the LT:

$$Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}+\frac{G\left(s\right)}{D\left(s\right)}$$

(14)

where $G\left(s\right)$ is the LT of the function $g\left(t\right)$. We denote the additional poles introduced by $G\left(s\right)$ as $s_v$ for $v=1,2,...,V$ (assuming that $s_k\ne s_v$, $\forall k,v$ and are of order $M=1$) with corresponding residues $c_v$. Applying the ILT and the residue theorem, we obtain the solution

$$y\left(t\right)=\sum _{k\in \mathbb{Z}}\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right)+\sum _{v=1}^V\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(\frac{c_v}{(s_v-a)e^{s_v\tau }-b}{e}^{s_{v}t}\right)\right),$$

(15)

where

$$\begin{array}{c}\hfill c_k=\frac{N\left(s_k\right)+G\left(s_k\right)}{D^{\prime }\left(s_k\right)},k\in \mathbb{Z}\end{array}$$

For more details about the LT, the interested readers are referred to [14]. In Section 4.1, we will explore the solutions of both cases (homogeneous and non-homogeneous).

3.2. Theory for NDDE Case

Next, we consider a simple NDDE case as shown below, where the solution is derived using the residue formula as described in Section 2.3:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =ay\left(t\right)+by(t-{\eta }_1)+c{y}^{\prime }(t-{\eta }_2)+g\left(t\right),\hfill & & t>0\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right),\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(16)

where $\tau =max({\eta }_1,{\eta }_2)$, and ${\eta }_1$, ${\eta }_2$$>0$. Considering the homogeneous case, and applying the LT to Equation (16) (using Equations (9) and (10)), we obtain the LT transform, $Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}$, of $y\left(t\right)$:

$$Y\left(s\right)=\frac{y\left(0\right)-cy(-{\eta }_2)+b{e}^{-s{\eta }_1}\underset{-{\eta }_1}{\overset{0}{\int }}y\left(v\right)e^{-sv}dv+cs{e}^{-s{\eta }_2}\underset{-{\eta }_2}{\overset{0}{\int }}y\left(v\right)e^{-sv}dv}{s-a-b{e}^{-s{\eta }_1}-cs{e}^{-s{\eta }_2}}.$$

(17)

Setting $D\left(s\right)=0$, we can find the poles $s_k$ of $Y\left(s\right)$. Under some general conditions, the poles $s_k$ are of order $M=1$. Then, applying the ILT, we obtain the solution

$$y\left(t\right)=\sum _{k\in \mathbb{Z}}\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right),$$

(18)

where the $c_k$s are the computed residues. We can analyze different cases depending on $g\left(t\right)$, and the signs and magnitudes of $a,b,c$ and ${\eta }_i$ for $i=1,2$. For further details, interested readers are referred to [14].

If the NDDE is non-homogeneous, similar properties hold as in the non-neutral case. In Section 4.2, we explore the solution behavior of different delay cases. Further details and analysis can be found in [14].

4. Comparison of MoS and LT Methods Using Symbolic Computation

Here, we present and discuss different examples of non-neutral DDEs and NDDEs introduced in [14] and apply the MoS and LT methodologies to each example to find the solution of such equations using the Maple and MATLAB software. First, the solutions obtained by each method will be compared between programs to highlight the variances in program computations. Second, the accuracy of the LT method will be evaluated against the analytical solution obtained by the MoS for each program. For this section, we introduce the following notation:

• $e_{MoS}\left(t\right)=|y_{Maple,MoS}\left(t\right)-y_{MATLAB,MoS}\left(t\right)|$: Absolute error in the solution $y\left(t\right)$ between Maple and MATLAB via the MoS.
• $e_{LT}\left(t\right)=|y_{Maple,LT}\left(t\right)-y_{MATLAB,LT}\left(t\right)|$: Absolute error in the solution $y\left(t\right)$ between Maple and MATLAB via the LT.
• $e_{Map}\left(t\right)=|y_{Maple,MoS}\left(t\right)-y_{Maple,LT}\left(t\right)|$: Absolute error in the solution $y\left(t\right)$ between the MoS and LT via Maple.
• $e_{MAT}\left(t\right)=|y_{MATLAB,MoS}\left(t\right)-y_{MATLAB,LT}\left(t\right)|$: Absolute error in the solution $y\left(t\right)$ between the MoS and LT via MATLAB.

4.1. Linear Non-Neutral DDE Examples

The two examples presented in this section demonstrate solutions formed from poles of order $M=1$. Figure 2 and Figure 3 and

Table 1. Computation time comparison for the non-neutral systems in Equations (19) and (20).

Example #	Program	Method	Time (s)
1	Maple	MoS	0.407
		LT	246.031
	MATLAB	MoS	1.335
		LT	64.914
2	Maple	MoS	0.485
		LT	269.828
	MATLAB	MoS	1.385
		LT	64.914

summarize the solution results, error comparisons between MoS and LT in each programming language, and the computation times for the corresponding example.

Example 1.

Let us consider the following DDE:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-15y\left(t\right)+15y(t-1),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=2-t-6t{(t+1)}^3,\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(19)

where $\tau =1$. Taking the LT of Equation (19), we obtain the quasi-polynomial $D\left(s\right)=(s+15)e^s-15$, which has a real pole at $s_0=0$ with residue $c_0=2.75$ (from Equation (4)). The remaining poles and residues are given by: $s_k=-15+W_k\left(15e^{15}\right)$, $k\in \mathbb{Z}\backslash \left\{0\right\}$, with property $s_{-k}=\overline{s_k}$, and $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$. Applying the ILT and Cauchy’s residue theorem, we obtain the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)=2.75+2\sum _{k=1}^{\infty }\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 1 via MoS (over $m=9$ intervals) and LT ($n=5000$) and in each software program is shown in Figure 2a. We notice that, as $t\to \infty$, $y\left(t\right)\to c_0=2.75$ since $a=-b$, $s_0=0$, and each Re$\left(s_k\right)<0$ for $k\in \mathbb{Z}\backslash \left\{0\right\}$. Figure 2b,c demonstrate the variances between the software when evaluating the solution via MoS (blue curve) or the LT (red curve). Figure 2d demonstrates the accuracy of the LT compared to the analytical MoS via Maple (blue curve) and MATLAB (red curve). Upon analyzing Figure 2b,d, we note that the MATLAB MoS solution significantly deviates from its LT solution and the Maple MoS solution. In the MoS process, a piecewise function is produced for each interval and it becomes harder to integrate and solve the sub-IVPs as the functions become more complicated with each step. In the LT solution, we see that, as t progresses, the variance fluctuations in the solution between programs decrease. Each plot shows that error spikes occur every $k\tau$. The reason is that, around these points, the sines of the truncated non-harmonic series do not play any significant role, while the cosine and exponential terms are essentially equal to 1. Therefore, the approximated series solution loses the ability to approximate the solution close to $t=k\tau$. The spikes can be reduced if we include more terms in the truncated series [14].

The computation times via each method and software are shown in Table 1. Given the number of terms to compute for the LT solution, the computation time was higher for both programs compared to the MoS solutions. However, with the use of MATLAB’s symbolic toolbox functions, the computation time for the LT solution was reduced by a factor of four compared to Maple for this particular example.

Example 2.

Let us consider the following DDE:

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-\frac{1}{5}y\left(t\right)-\frac{5}{8}y(t-3),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)={(t+3/2)}^2,\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(20)

where $\tau =3$. Taking the LT of Equation (20), we obtained the quasi-polynomial $D\left(s\right)=(s-1/5)e^{3s}-5/8$, which has no real poles since $D\left(\ln (-b\tau )/\tau \right)>0$ (see [14]). The complex poles are given by: $s_k=-\frac{1}{5}+\frac{1}{3}{W}_k\left(-\frac{15}{8}{e}^{3/5}\right)$, $k\in \mathbb{Z}$, with property $s_{-(k+1)}=\overline{s_k}$, and corresponding residues $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$. Applying the ILT and Cauchy’s residue theorem, we obtained the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)=2\sum _{k=1}^{\infty }Re\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 2 via MoS (over $m=9$ intervals) and LT ($n=5000$) and in each software program is shown in Figure 3a. We noticed that, as $t\to \infty$, $y\left(t\right)\to 0$ since Re$\left(s_k\right)<0$$\forall k\in \mathbb{Z}$. Figure 3b,c demonstrates the variances between the software when evaluating the solution via MoS (blue curve) or the LT (red curve). Figure 3d demonstrates the accuracy of the LT compared to the analytical MoS via Maple (blue curve) and MATLAB (red curve).

Compared to Example 1, Figure 3b,d demonstrates that the performance error in MATLAB’s MoS solution, relative to the Maple MoS solution and the MATLAB LT solution, decreased. One notable difference between Examples 1 and 2 is the reduced complexity in the history function. In Example 1, the history function was a fifth-order polynomial, whereas, here, the history is a second-order polynomial. In both examples, the LT solution is truncated to 5000 terms. Note that truncation creates ringing in the solution near the history function endpoints. In Figure 3d, the LT error decreases exponentially for both software, while, in Figure 2d, we notice that the error has a maximum peak at around $t=8$. In some cases, increasing the number of terms reduces the error further and improves the history function approximation, but at the cost of computational complexity and time.

The computation times via each method and software are shown in Table 1. We noticed that Examples 1 and 2 exhibit similar times since we used the same solution structure in both programs, compatible functions between programs, number of terms to be computed, and number of intervals to integrate over. Given the number of terms to be computed for the LT solution, the computation time was higher for both programs compared to the MoS solutions. However, with the use of MATLAB’s symbolic toolbox functions, the computation time for the LT solution was reduced by 4x compared to Maple for both examples.

4.2. Linear NDDE Examples

The NDDE examples to be discussed conform to the general form in Equation (16). These five examples demonstrate solutions formed from poles of order $M=1$. Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 and

Table 2. Computation time comparison for the NDDE systems in Equations (21)–(25).

Example #	Program	Method	Time (s)
3	Maple	MoS	0.203
		LT	84.219
	MATLAB	MoS	0.838
		LT	870.978
4	Maple	MoS	0.296
		LT	15.470
	MATLAB	MoS	2.022
		LT	32.663
5	Maple	MoS	0.453
		LT	17.110
	MATLAB	MoS	1.572
		LT	29.052
6	Maple	MoS	0.438
		LT	42.344
	MATLAB	MoS	1.130
		LT	284.252
7	Maple	MoS	0.297
		LT	89.093
	MATLAB	MoS	1.280
		LT	349.196

summarize the solution results, error comparisons between MoS and LT in each programming language, and the computation times for the corresponding example.

Example 3.

Consider

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =y\left(t\right)+y(t-1)-\frac{1}{4}{y}^{\prime }(t-1),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=-t,\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(21)

where ${\eta }_1=1$, ${\eta }_2=1$, and $\tau =max\{{\eta }_1,{\eta }_2\}=1$. This example was utilized in [36] to evaluate the performance of a numerical method for solving linear NDDEs using the generalized Lambert W function [55,56]. Applying MoS for a few partitions, we obtained the following analytical solution:

$$\begin{array}{c}\hfill y\left(t\right)=\left\{\begin{array}{cc}-t,\hfill & [-\tau ,0],\hfill \\ t+\frac{1}{4}{e}^t-\frac{1}{4},\hfill & [0,1],\hfill \\ -t+\frac{1}{2}+\left(\frac{3}{16}t+\frac{17}{16}+\frac{1}{4}e\right)e^{t-1},\hfill & [1,2],\hfill \\ \cdots \hfill & \cdots \hfill \end{array}\right.\end{array}$$

Taking the LT of Equation (21), we obtained the quasi-polynomial $D\left(s\right)=s-1-e^{-s}+\frac{1}{4}s{e}^{-s}$, which has a real pole at $s_0\approx 1.2084$ with residue $c_0\approx 0.3837$ (from Equation (4)). The remaining poles and residues are given by the sequence ${\left\{s_k\right\}}_{k=-\infty }^{\infty }$ for $k\in \mathbb{Z}\backslash \left\{0\right\}$$(s_{-k}=\overline{s_k})$ and $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$, respectively. For large $j\ge k$, these poles approach $s_j=\frac{\ln \left(\right|b\left|\right)+(2j+1)\pi i}{\tau }\approx -1.3863+(2j+1)\pi i$ for $j\in \mathbb{N}$ (see [14]). Since $s_0>0$, then, from the stability analysis, the solution will be unstable. Applying the ILT and Cauchy’s residue theorem, we obtained the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)\approx 0.3837e^{1.2084t}+2\sum _{k=1}^{\infty }Re\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 3 via MoS (over $m=5$ intervals) and LT ($n=2000$) and in each software program is shown in Figure 4a. Figure 4b,c shows that, as $t\to \infty$, the variance between each software increases exponentially. Figure 4d indicates that the accuracy of the LT solution improves exponentially as t increases for both programs, and with noticeable spikes every $k\tau$. Compared to the non-neutral examples, NDDE solutions have a higher complexity and, thus, one software may be more beneficial to use vs. another.

Example 4.

Consider

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-2y\left(t\right)+y(t-2)+\frac{1}{2}{y}^{\prime }(t-2),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=\sin \left(\pi t\right),\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(22)

where ${\eta }_1=2$, ${\eta }_2=2$, and $\tau =max\{{\eta }_1,{\eta }_2\}=2$. This second example was first presented in [57] to study the superconvergence of continuous Galerkin finite element methods for linear NDDEs. Applying the MoS, the following analytical solution of $y\left(t\right)$ over the interval $[-\tau ,m\tau ]$ was obtained, where m is the number of intervals:

$$\begin{array}{c}\hfill y\left(t\right)=\left\{\begin{array}{cc}\sin \left(\pi t\right)\hfill & [-\tau ,0]\hfill \\ \sin \left(\pi t\right)/2\hfill & [0,2]\hfill \\ \sin \left(\pi t\right)/4\hfill & [2,4]\hfill \\ ⋮\hfill & ⋮\hfill \\ \sin \left(\pi t\right)/\left(2^m\right)\hfill & \left[\right(m-1)\tau ,m\tau ].\hfill \end{array}\right.\end{array}$$

The characteristic equation of the system is given by $D\left(s\right)=s+2-\frac{1}{2}s{e}^{-2s}-e^{-2s}$, which has a real pole at $s_0\approx -0.3466$ with residue $c_0\approx -0.0786$. The remaining poles and residues are given by the sequence ${\left\{s_k\right\}}_{k=-\infty }^{\infty }$ for $k\in \mathbb{Z}\backslash \left\{0\right\}$$(s_{-k}=\overline{s_k})$ and $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$, respectively. These poles are exactly given by $s_k=\frac{\ln \left(c\right)+2k\pi i}{\tau }\approx -0.3466+k\pi i$ for $k\in \mathbb{N}$ (see [14]). Since $s_0<0$ and Re$\left(s_k\right)<0$$\forall k$, then, from the stability analysis, the solution will be stable. Applying the ILT and Cauchy’s residue theorem, we obtained the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)\approx -0.0786e^{-0.3466t}+2\sum _{k=1}^{\infty }Re\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 4 via MoS (over $m=6$ intervals) and LT ($n=500$) and in each software program is shown in Figure 5a. Figure 5b,c shows that the variance between each software is of order ${10}^{-16}$ and, as $t\to \infty$, the error decreases with peak errors every $t=k$ for $k\in \mathbb{N}$. Figure 5d indicates that the accuracy of the LT solution improves as t increases for both programs, and has an interesting spike located around $\pi \tau$.

Example 5.

Let us consider

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-\frac{3}{2}y\left(t\right)-\frac{9}{10}{y}^{\prime }(t-2)+6,\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=-\frac{5}{6}(t+2)t^3,\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(23)

where ${\eta }_2=2$ and $\tau ={\eta }_2$. The third example does not contain a state-dependent delay term. This equation is governed by the characteristic quasi-polynomial given by $D\left(s\right)=s+\frac{3}{2}+\frac{9}{10}s{e}^{-2s}$, which has a real pole at $s_0\approx -0.4602$ with residue $c_0\approx -2.2994$. The other poles are given by the sequence ${\left\{s_k\right\}}_{k=1}^{\infty }$$(s_{-k}=\overline{s_k})$. For large $j\ge k$, these poles approach $s_j=\frac{\ln \left(\right|c\left|\right)+(2j+1)\pi i}{\tau }\approx -0.0527+(2j+1)\pi i$ for $j\in \mathbb{N}$ (see [14]). The corresponding residues are given by $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$. Applying the ILT and Cauchy’s residue theorem, we obtained the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)\approx -2.2994e^{-0.4602t}+2\sum _{k=1}^{\infty }Re\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 5 via MoS (over $m=10$ intervals) and LT ($n=500$) and in each software program is shown in Figure 6a. Figure 6b,c shows that the variance between each software is of order ${10}^{-5}$ via MoS and order ${10}^{-15}$ via LT. From Figure 6d, as $t\to \infty$, the error decreases with peak errors every $t=k$ for $k\in \mathbb{N}$. This indicates that the accuracy of the LT solution improves as t increases for both programs, and has an interesting spike located around $\pi \tau$.

Example 6.

Let us consider

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-6y\left(t\right)-6y(t-1)+\frac{9}{10}{y}^{\prime }(t-2),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=-t(1+t),\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(24)

where ${\eta }_1=1$, ${\eta }_2=2$, and $\tau =max\{{\eta }_1,{\eta }_2\}={\eta }_2$. The fourth example has the characteristic quasi-polynomial given by $D\left(s\right)=s+6-\frac{9}{10}s{e}^{-2s}+6e^{-s}$, which has no real poles. The complex poles are given by the sequence ${\left\{s_k\right\}}_{k=1}^{\infty }$$(s_{-k}=\overline{s_k})$. For large $j\ge k$, these poles approach $s_j=\frac{\ln \left(c\right)+2j\pi i}{\tau }\approx -0.0527+j\pi i$ for $j\in \mathbb{N}$ (see [14]). The corresponding residues are given by $c_k=\frac{N\left(s_k\right)}{D^{\prime }\left(s_k\right)}$. Applying the ILT and Cauchy’s residue theorem, we obtained the LT solution of $y\left(t\right)$:

$$\begin{array}{c}y\left(t\right)\approx 2\sum _{k=1}^{\infty }Re\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right).\end{array}$$

The solution of Example 6 via MoS (over $m=7\tau$ intervals) and LT ($n=500$) and in each software program is shown in Figure 7a. Figure 7b,c shows that the variance between each software is of order ${10}^{-9}$ via MoS and order ${10}^{-16}$ via LT. We see similar deviances in the MoS solutions compared to Example 5 between each software. A larger timespan would be needed to deduce conclusions on the similarities of the LT solution between each software. From Figure 6d, as $t\to \infty$, the error decreases, with peak errors every $t=k\tau$ for $k\in \mathbb{N}$.

Example 7.

Let us consider

$$\begin{array}{cccc}\hfill y^{\prime }\left(t\right)& =-6y\left(t\right)+6y(t-2)+\frac{2}{3}{y}^{\prime }(t-1),\hfill & & t>0,\hfill \\ \hfill y\left(t\right)& ={\varphi }_0\left(t\right)=4+t,\hfill & & t\in [-\tau ,0],\hfill \end{array}$$

(25)

where ${\eta }_1=2$, ${\eta }_2=1$, and $\tau =max\{{\eta }_1,{\eta }_2\}={\eta }_1$. The fifth example in Equation (25) has the characteristic quasi-polynomial given by $D\left(s\right)=s+6-\frac{2}{3}s{e}^{-s}-6e^{-2s}$, which has a real pole at $s_0=0$ with residue $c_0\approx 3.0811$. The remaining poles are divided into two sequences (see [14]). One sequence of poles, ${\left\{s_k\right\}}_{k=1}^{\infty }$$(s_{-k}=\overline{s_k})$, lies relatively close to the imaginary axis, which, for large $j\ge k$, approaches $s_j=\frac{\ln \left(c\right)+2j\pi i}{{\eta }_2}\approx -0.4055+2j\pi i$ for $j\in \mathbb{N}$. The second sequence of poles, ${\left\{s_{k^{\prime }}\right\}}_{k^{\prime }=1}^{\infty }$$(s_{-k^{\prime }}=\overline{s_{k^{\prime }}})$, is offset from the imaginary axis, which, for large $j^{\prime }\ge k^{\prime }$, approaches $s_{j^{\prime }}=\frac{W_{j^{\prime }}\left(\frac{b}{c}\left({\eta }_2-{\eta }_1\right)\right)}{{\eta }_1-{\eta }_2}$ for $j^{\prime }\in \mathbb{N}$. The corresponding residues are given by $c_{k,k^{\prime }}=\underset{s\to s_k,s_{k^{\prime }}}{lim}\left(\frac{N\left(s\right)}{1+\frac{2}{3}(s-1)e^{-s}+12e^{-2s}}\right)$. Then, the solution is of the form

$$\begin{array}{c}y\left(t\right)\approx 3.0811+2\sum _{k=1}^{\infty }\mathrm{Re}\left(\underset{s\to s_k}{lim}\left(c_k{e}^{s_{k}t}\right)\right)+2\sum _{k^{\prime }=1}^{\infty }Re\left(\underset{s\to s_{k^{\prime }}}{lim}\left(c_{k^{\prime }}{e}^{s_{k^{\prime }}t}\right)\right).\end{array}$$

The solution of Example 7 via MoS (over $m=7\tau$ intervals) and LT ($n=500$) and in each software program is shown in Figure 8a. Figure 8b,c shows that the variance between each software is of order ${10}^{-9}$ via MoS and order ${10}^{-15}$ via LT. We see similar deviances in the MoS solutions near the boundary intervals compared to Example 5. Relative to previous NDDE examples, Example 7 has a larger error in the LT solution vs. the analytical MoS via both software, with an order of magnitude ${10}^{-3}$ as shown in Figure 8d, despite the simplicity of the history function. Nevertheless, we can observe that, as $t\to \infty$, the error decreases.

The computation times via each method and software are shown in Table 2. We make the following notes regarding the time computation results:

• From Maple MoS solutions, the order of time computation is as follows:
  Example 5 > Example 6 > Example 7 > Example 4 > Example 3.
• From Maple LT solutions, the order of time computation is as follows:
  Example 7 > Example 3 > Example 6 > Example 5 > Example 4.
• From MATLAB MoS solutions, the order of time computation is as follows:
  Example 4 > Example 5 > Example 7 > Example 6 > Example 3.
• From MATLAB LT solutions, the order of time computation is as follows:
  Example 3 > Example 7 > Example 6 > Example 4 > Example 5.

From these results, we notice that the computation time is dependent on the complexity of the history function, the number of terms used in the LT solution, the number of intervals used in the MoS solution, and the magnitude of the parameters of the NDDE. Regarding MoS solutions, Examples 4 and 5 both have a complicated history function (one is a fourth-order polynomial and the other is periodic), and we observe that both programs took longer to compute the solutions. Regarding LT solutions, Examples 3 and 7 took the most time to compute via both software, even though Example 1 had the most terms. The differences between these examples is that Example 7 consisted of an extra sequence of poles due to the different delays in the state-dependent term and the derivative-of-state-dependent term. Nevertheless, Example 3 was found to be the quickest to compute via MoS for both softwares, whereas Examples 4 and 5 were found to be the quickest to compute via LT. We also note that, in all cases, Maple’s computation time was found to significantly outrank MATLAB’s computation time.

We notice that Examples 3 and 4 exhibit similar times since we used the same solution structure in both programs, compatible functions between programs, number of terms to be computed, and number of intervals to integrate over. Given the number of terms to be computed for the LT solution, the computation time was higher for both programs compared to the MoS solutions. However, with the use of MATLAB’s symbolic toolbox functions, the computation time for the LT solution was reduced by 4x compared to Maple for both examples.

5. Conclusions

In this paper, we investigated the symbolic mathematical computational performance of the software program Maple and the programming language MATLAB when applying the method of steps (MoS) and the Laplace transform (LT) methods to solve neutral and retarded linear DDEs. Of particular interest was comparing how these programs performed. The exact, analytic MoS solution was used to aid in this comparison. However, the main reason for why most of the graphs only include the first 5–10 time-intervals is due to the fact that the method of steps often could not compute the solution beyond these intervals. The Laplace solution is an infinite series, which requires one to compute the relevant poles along with their associated residues. We computed the analytical solutions using Maple and MATLAB in order to take advantage of their respective symbolic computation features. The accuracy of the Laplace method solutions were measured by comparing them with those obtained by the method of steps. The Laplace method produces an analytical solution that can be evaluated at any particular value of time. The results obtained emphasize that symbolic computational software is a powerful tool, capable of computing analytic solutions for linear delay differential equations, including neutral ones. From a computational viewpoint, we found that the computation time is dependent on the complexity of the history function, the number of terms used in the LT solution, the number of intervals used in the MoS solution, and the DDE parameters. We found that, for linear non-neutral DDEs, MATLAB symbolic computations were faster than the ones from Maple. However, for the linear neutral DDEs, which are often more difficult to solve, Maple was faster. Regarding accuracy, it is important to point out that we compared Maple and MATLAB via MoS and LT. The differences between both software were smaller for the LT in comparison to the MoS. In addition, the accuracy of the LT solutions using Maple were, in a few cases, slightly better than the ones generated by MATLAB. However, both were reliable, but we can mention that Maple lightly surpasses MATLAB in symbolic computation due to its symbolic capabilities/nature. Finally, we would like to mention some limitations of the LT irrespective of the software choice. The LT method is only applicable for linear DDEs due to the Laplace transform definition. Thus, nonlinear DDEs should be solved by other means. In addition, we implemented the LT for scalar linear DDEs with constant coefficients and a finite number of delays. However, the LT method can also be extended to solve first-order linear DDE systems. In this work, we applied the LT method to homogeneous and non-homogenous DDEs. Future works that are in progress include solving linear non-homogeneous DDEs with Dirac delta function inputs, with one related goal being to develop a Green’s function type kernel.