 Research
 Open Access
 Published:
Numerical algorithm for nonlinear delayed differential systems of nth order
Advances in Difference Equations volume 2019, Article number: 26 (2019)
Abstract
The purpose of this paper is to propose a semianalytical technique convenient for numerical approximation of solutions of the initial value problem for pdimensional delayed and neutral differential systems with constant, proportional and time varying delays. The algorithm is based on combination of the method of steps and the differential transformation. Convergence analysis of the presented method is given as well. Applicability of the presented approach is demonstrated in two examples. A system of pantograph type differential equations and a system of neutral functional differential equations with three types of delays are considered. The accuracy of the results is compared to those obtained by the Laplace decomposition algorithm, the residual power series method and Matlab package DDENSD. A comparison of computing time is presented, too, showing reliability and efficiency of the proposed technique.
Introduction
Systems of functional differential equations (FDEs), in particular delayed or neutral differential equations, are often used to model processes in the real world. To give some examples, we mention models in population dynamics [1], neuromechanics [2], machine tool vibrations [3], etc. Further models and details can be found, for instance, in monographs [4] and [5].
Semianalytical methods expressing solutions to problems with delays in a series form have been studied in the last two decades. Methods such as the variational iteration method (VIM) [6], Adomian decomposition method (ADM) [7], homotopy perturbation method (HPM) [8], homotopy analysis method (HAM) [9] and also methods based on the Taylor theorem such as the differential transformation (DT) [10], Taylor collocation method [11] and Taylor polynomial method [12] have been developed to approximate solutions to different problems for FDEs. Other ways to use the series approach in solving FDEs are, e.g., the method of polynomial quasisolutions [13, 14], finite difference methods [15, 16], and the functional analytic technique (FAT) [17, 18].
The main aim of the work is to apply a combination of the method of steps and DT as a convenient tool for finding an approximate solution to the initial value problem for functional differential systems used in dynamical models. Convergence analysis and error estimates of the method are investigated as well. We give some experimental results in Sect. 4 to show that the algorithm produces reliable results with the same or better efficiency than the reference methods.
Methods
The main idea of our approach is to combine the differential transformation and general method of steps.
The differential transformation has been, and still is, an active research topic during the last years. As examples of recently published results, we mention research papers [19,20,21,22,23]. These papers among other publications contain new algorithms and their applications to solving different problems involving differential equations.
Definition 1
The differential transformation of a real function \(u(t)\) at a point \(t_{0} \in \mathbb{R}\) is \(\mathcal{D} \{ u(t) \} [t_{0}] = \{ U(k)[t _{0}] \}_{k=0}^{\infty }\), where the kth component \(U(k)[t_{0}]\) of the differential transformation of the function \(u(t)\) at \(t_{0}\) is defined as
assuming that the original function \(u(t)\) is analytic.
Definition 2
The inverse differential transformation of \(\{ U(k)[t_{0}] \}_{k=0} ^{\infty }\) at \(t_{0}\) is defined as
In applications, the function \(u(t)\) is usually expressed in the form of finite series
In Sect. 4, we use the following transformation formulas, which are derived from definitions (1), (2) and proved in [24].
Lemma 1
Assume that \(W(k)\), \(U(k)\) and \(U_{i} (k)\) are the kth components of the differential transformations of functions \(w(t)\), \(u(t)\) and \(u_{i} (t)\), \(i=1,2\), at \(t_{0} \in \mathbb{R}\), respectively, and let \(q, q_{j} \in (0,1)\), \(j=1,2\). Moreover, assume that \(t_{0}=0\). Denote \(\mathbb{N}_{0} = \mathbb{N} \cup \{ 0 \}\).
Remark 1
Transformation formulas for shifted arguments \(w(t)=u(ta)\) are often proved and applied in papers. However, using these formulas when solving initial value problems for delayed differential equations is not convenient since the uniqueness of solutions is violated. The reason is that the values of the initial vector function for \(t < 0\) are not taken into account.
One of the drawbacks of the common approach to the differential transformation is that there is no use of direct transformation formulas for equations with nonlinear terms containing unknown function \(u(t)\), for instance, \(f(u)= \operatorname {e}^{\cos {u}}\) or \(f(u) = \sqrt{1+u^{4}}\).
Fortunately, the corresponding transformations can be calculated using the Adomian polynomials \(A_{n}\), in which each solution \(u_{i}\) is replaced by the corresponding components \(U_{i} (k)\) of the differential transformation \(\{ U_{i} (k) \}_{k=0}^{\infty }\), see [25]. Suppose that \(F(k)\) is the kth component of the differential transformation of a nonlinear term \(f(u)\), then
Recently, it turned out that there is another way to work with nonlinearities in DT [26].
The second method, namely the method of steps, enables us to replace the terms involving constant or timedependent delays by the initial vector function and its derivatives. Then the original initial value problem for a system of delayed or neutral differential equations is simplified to the initial problem for a system of ordinary differential equations. Details on the method of steps can be found, e.g., in monographs [4, 5, 27].
Results
The subject of our interest is a system of p functional differential equations of nth order with multiple delays \(\alpha _{1} (t),\ldots, \alpha _{r} (t)\) in the following form:
where \(\mathbf{u}^{(n)} (t) =(u_{1}^{(n)}(t),\ldots, u_{p}^{(n)}(t))^{T}\), \({\mathbf{u}^{(k)} (t) = (u_{1}^{(k)}(t),\ldots, u_{p}^{(k)}(t))^{T}}\), \(k=0,1,\ldots, n1\) and \(\mathbf{f} = (f_{1},\ldots, f_{p})^{T}\) are pdimensional vector functions, \(\mathbf{u} _{i}(\alpha _{i}(t))= (\mathbf{u}(\alpha _{i}(t)),\mathbf{u}'(\alpha _{i}(t)),\ldots, \mathbf{u}^{(m_{i})}(\alpha _{i}(t)))\) are \((m_{i} \cdot p)\)dimensional vector functions, \(m_{i} \leq n\), \(i=1,2,\ldots,r\), \(r \in \mathbb{N}\) and \(f_{j} \colon [0,\infty ) \times \mathbb{R}^{np} \times \mathbb{R}^{\omega p}\) are continuous real functions for \(j=1,2,\ldots,p\), where \(\omega = \sum_{i=1}^{r} m_{i}\).
We consider three types of delays \(\alpha _{i}\):

1.
\(\alpha _{i}(t) = q_{i} t\), where \(q_{i} \in (0,1)\) (proportional delay).

2.
\(\alpha _{i}(t) = t\tau _{i}\), where \(\tau _{i}>0\) is a real constant (constant delay).

3.
\(\alpha _{i}(t)=t\tau _{i}(t)\), where \(\tau _{i}(t) \geq \tau _{i0}>0\) for \(t>0\) is a real function (timedependent or timevarying delay).
Let \(t^{*}= \min_{1\leq i \leq r} \{\inf_{t>0} (\alpha _{i}(t) ) \} \leq 0\), \(m= \max \{m_{1},m_{2},\ldots,m_{r}\} \leq n\). In the case \(m=n\), system (5) is a neutral system, otherwise it is a delayed differential system.
If \(t^{*}<0\), an initial vector function \(\varPhi (t) = (\phi _{1}(t),\ldots,\phi _{p}(t))^{T}\) must be assigned to system (5) on the interval \([t^{*}, 0]\). Moreover, we assume that \(\phi _{j}(t) \in C ^{n}([t^{*},0],\mathbb{R})\) for \(j=1,\ldots,p\).
We look for a solution of system (5) with the following initial conditions:
and the initial vector function \(\varPhi (t)\) on interval \([t^{*}, 0]\) satisfying
We solve initial value problem (5), (6) and (7) subject to the following hypotheses:

(H1)
The functions \(f_{j}\), \(j=1,\ldots, p\) are analytic in \([0,T^{*}] \times \mathbb{R}^{np} \times \mathbb{R}^{\omega p}\).

(H2)
The initial value problem (5), (6) and (7) has a unique solution on some interval \([0, T^{*}]\).
Remark 2
Hypothesis (H2) is valid, for example, if the delay functions \(\alpha _{i}\) are Lipschitz continuous on \([0,T^{*}]\), the functions \(\phi _{j}, \phi '_{j},\ldots, \phi _{j}^{(n)}\) are Lipschitz continuous on \([t^{*},0]\), and the functions \(f_{j}\) are continuous with respect to t on \([0,T^{*}]\) and Lipschitz continuous with respect to the rest of the variables on \(\mathbb{R}^{np} \times \mathbb{R}^{\omega p}\). More details and other types of sufficient conditions for existence of a unique solution can be found in [5, Sects. 3.2 and 3.3], or [27, Sect. 2.2].
We start with the method of steps. We substitute the initial vector function \(\varPhi (t)\) and its derivatives in all places where the unknown functions with constant or timedependent delays and derivatives of those functions take place. This turns the delayed system (5) into a system of ordinary differential equations or differential equations with proportional delays in the case when system (5) contains proportional delays.
For example, if \(\alpha _{1} (t) = t\tau _{1}\), \(\alpha _{2} (t) = t \tau _{2}\), \(\alpha _{3} (t) = q_{3} t\) and \(\alpha _{4} (t) = t\tau _{4}(t)\), applying the method of steps changes (5) into the system
where
and \(m_{l} \leq n\) for \(l=1,2,3,4\). Then we transform the initial conditions (6). Definition (1) gives
After applying the differential transformation, the initial value problem for a system of FDEs is reduced to a system of recurrence algebraic relations
Solving this recurrence and then using the inverse transformation (2), we get an approximate solution of system (5) in the series form
If \(t^{*}<0\), we denote \(t_{\alpha _{i}} = \inf \{ t: \alpha _{i} (t)>0 \}\) and \(t_{\alpha } =\min_{1 \leq i \leq r} \{ t_{\alpha _{i}}: t_{\alpha _{i}} \neq 0 \}\). Then the approximate solution \(\mathbf{u}(t)\) is valid on the intersection of its convergence interval and the interval \([0,T^{*}] \cap [0, t_{\alpha }]\), whereas \(\mathbf{u}(t) = \varPhi (t)\) on the interval \([t^{*},0]\). If \(t^{*}=0\), the approximate solution \(\mathbf{u}(t)\) is valid on the intersection of its convergence interval with \([0,T^{*}]\).
Now we formulate and prove two theorems on convergence and an error estimate of the approximate solution to the studied problem obtained using the differential transformation.
Theorem 1
Let hypotheses (H1) and (H2) be valid and denote \(\mathbf{F}_{k}(t) = \mathbf{U}(k)t^{k}\). If there exist a constant δ, \(0<\delta < 1\), and \(k_{0} \in \mathbb{N}\) such that \(\Vert \mathbf{F}_{k+1}(t)\Vert \leq \delta \Vert \mathbf{F}_{k}(t)\Vert \) for all \(k \geq k_{0}\), then the series \(\sum_{k=0}^{\infty } \mathbf{F}_{k}(t)\) converges to a unique solution on the interval \(J=[0,\gamma ]\), \(\gamma \leq T^{*}\).
Proof
Denote \(C^{n}(J)\) the Banach space of vectorvalued functions \(\mathbf{h}(t) = (h_{1}(t),h_{2}(t),\ldots, h_{p}(t))^{T}\) with continuous derivatives up to order n and norm
Denote
Now it is sufficient to prove that sequence \(\{\mathbf{S}_{l} \} \) is a Cauchy sequence in the Banach space \(C^{n}(J)\). Due to
for every \(l,m \in \mathbb{N}\), \(l \geq m >n_{0}\), we get
Since \(0<\delta <1\), it follows that
Therefore, \(\{\mathbf{S}_{l} \} \) is a Cauchy sequence in the Banach space \(C^{n}(J)\) and the proof is complete. □
Theorem 2
Suppose that the assumptions of Theorem 1 are valid. Then for the truncated series \(\sum_{k=0}^{m}\mathbf{F}_{k}(t)\) the following error estimate holds:
for any \(m_{0} \geq 0\), \(m \geq m_{0}\).
Proof
Without loss of generality, we can choose \(m_{0} \geq n\), where n is the order of system (5). From inequality (10) we have
for \(l \geq m \geq m_{0}\). From \(0 <\delta <1\) it follows \((1 \delta ^{lm})< 1\). Hence inequality (11) can be reduced to
Here we use the fact that for \(l \rightarrow \infty \), \(\mathbf{S} _{l} \rightarrow \mathbf{u}(t)\), and so the proof is complete. □
Remark 3
Recent results on error estimates and convergence of Taylor series can be found, e.g., in [28].
Applications and discussion
As the first application, we have chosen the initial value problem, which has been solved in [29] using the Laplace decomposition method (LDM) and in [30] using the residual power series method (RPSM).
Example 1
We are looking for a solution of a 3dimensional system of pantograph equations
subject to the initial conditions
Since system (12) contains proportional delays only, we do not have to use the method of steps. Applying DT formulas in Lemma 1 to (12), we get a system of recurrence relations
From the initial conditions we have \(U_{1}(0)=1\), \(U_{2}(0)=0\), \(U_{3}(0)=0\). Solving system (14), we get
For \(k \geq 2\), we find
Application of the inverse differential transformation (2) gives a solution to (12), (13) in the form
When \(N \rightarrow \infty \), the series converge to the Taylor expansions of the closedform solutions
Comparison of absolute errors of the presented DT technique with LDM and RPSM for \(N=2\) is done in Tables 1, 2 and 3. We see that DT and RPSM produce the same results, which are close to the values of the closed form solutions, whereas LDM shows significant deviations. We obtain similar results when comparing computing times, see Tables 4, 5 and 6.
Remark 4
In [29], the authors used LDM and obtained only approximate solutions of the initial value problem (12), (13). Applying RPSM, the authors were able to find closedform solutions in [30]. However, the calculations are too complicated, and the residual functions (RPSM) and initial guesses (LDM) contain analytical forms of functions sin and cos, which means that these methods are not convenient for use in a purely numerical software.
As the second application, we have chosen a system with all three types of delays considered to show reliability and efficiency of the proposed approach in solving difficult tasks.
Example 2
Let us solve a nonlinear system of neutral delayed differential equations
with initial functions
for \(t \in [2,0]\), and initial conditions
Following the method of steps, we get
System (18) cannot be solved by the classical DT approach because of the nonlinear term \(f(u) = \sqrt[3]{u_{1}^{2}}\), hence we apply modified Adomian formula for differential transformation components to the nonlinear term \(f(u)\). Applying DT to (18), we get the system
where \(F_{1}(k)\) is the kth component of the transformed function \(f(u) = \sqrt[3]{u_{1}^{2}}\). Applying formula (4) and the transformed initial conditions
we obtain
Solving the system of recurrence relations (19)–(20), we get
Applying the inverse differential transformation, we obtain an approximate solution to the initial value problem (15)–(17):
As we do not know the exact solution of the given problem, we are limited to comparing approximate solutions. Comparison of values obtained by the proposed approach and values obtained by Matlab package DDENSD in Table 7 shows good correspondence between the results. Comparing computing times in Table 8, we can see that the presented method produces reliable results much faster than Matlab package DDENSD.
Remark 5
System (15) contains all three types of delay which were considered in this paper. Moreover, it contains a term which is nonlinear (nonpolynomial) in the dependent variable \(u_{1}\). In this sense, the present paper contains more complicated systems in applications than papers about other semianalytical methods like VIM [6], ADM [7] or HPM [8].
Conclusions
The approach presented in this paper is an effective semianalytical technique, convenient for numerical approximation of a unique solution to the initial value problem for systems of functional differential equations, in particular delayed and neutral differential equations. Considering systems of equations with three types of delays brings a generalization with respect to the problem studied in [20]. The comparison of results was done against the Laplace decomposition method, residual power series method and Matlab package DDENSD. The need of computational work is reduced compared to the other methods. The differential transformation algorithm gives an approximate solution which is in good concordance with reference results produced by Matlab. Under certain circumstances, it is possible to identify the unique solution to the initial value problem in closed form. Further steps can be done in the development of the presented technique for systems with distributed and state dependent delays.
References
 1.
Györi, I.: Oscillation and comparison results in neutral differential equations and their applications to the delay logistic equation. Comput. Math. Appl. 18(10–11), 893–906 (1989)
 2.
Insperger, T., Milton, J., Stepan, G.: Semidiscretization for time delayed neural balance control. SIAM J. Appl. Dyn. Syst. 14(3), 1258–1277 (2015)
 3.
KalmarNagy, T., Stepan, G., Moon, F.C.: Subcritical Hopf bifucration in the delay equation model for machine tool vibrations. Nonlinear Dyn. 26, 121–142 (2011)
 4.
Hale, J.K., Verduyn Lunel, S.M.: Introduction to Functional Differential Equations. Springer, New York (1993)
 5.
Kolmanovskii, V., Myshkis, A.: Introduction to the Theory and Applications of Functional Differential Equations. Kluwer Academic, Dordrecht (1999)
 6.
Chen, X., Wang, L.: The variational iteration method for solving a neutral functionaldifferential equation with proportional delays. Comput. Math. Appl. 59, 2696–2702 (2010)
 7.
BlancoCocom, L., Estrella, A.G., AvilaVales, E.: Solving delay differential systems with history functions by the Adomian decomposition method. Appl. Math. Comput. 218, 5994–6011 (2013)
 8.
Shakeri, F., Dehghan, M.: Solution of delay differential equations via a homotopy perturbation method. Math. Comput. Model. 48, 486–498 (2008)
 9.
Duarte, J., Januario, C., Martins, N.: Analytical solutions of an economic model by the homotopy analysis method. Appl. Math. Sci. 10(49), 2483–2490 (2016)
 10.
Rebenda, J., Šmarda, Z.: A semianalytical approach for solving nonlinear systems of functional differential equations with delay. In: Simos, T.E. (ed.) 14th International Conference of Numerical Analysis and Applied Mathematics (ICNAAM 2016). AIP Conference Proceedings, vol. 1863, p. 530003. AIP Publishing, Melville (2017)
 11.
Bellour, A., Bousselsal, M.: Numerical solution of delay integrodifferential equations by using Taylor collocation method. Math. Methods Appl. Sci. 37, 1491–1506 (2014)
 12.
Sezer, M., AkyuzDascioglu, A.: Taylor polynomial solutions of general linear differentialdifference equations with variable coefficients. Appl. Math. Comput. 174, 753–765 (2006)
 13.
Cherepennikov, V.B., Ermolaeva, P.G.: Smooth solutions of an initialvalue problem for some differential difference equations. Numer. Anal. Appl. 3, 174–185 (2010)
 14.
Cherepennikov, V.B.: Numerical analytical method of studying some linear functional differential equations. Numer. Anal. Appl. 6, 236–246 (2013)
 15.
Jain, R.K., Agarwal, R.P.: Finite difference method for second order functional differential equations. J. Math. Phys. Sci. 7(3), 301–3016 (1973)
 16.
Agarwal, R.P., Chow, Y.M.: Finitedifference methods for boundaryvalue problems of differential equations with deviating arguments. Comput. Math. Appl. 12A(11), 1143–1153 (1986)
 17.
Petropoulou, E.N., Siafarikas, P.D., Tzirtzilakis, E.E.: A “discretization” technique for the solution of ODEs. J. Math. Anal. Appl. 331, 279–296 (2007)
 18.
Petropoulou, E.N., Siafarikas, P.D., Tzirtzilakis, E.E.: A “discretization” technique for the solution of ODEs II. Numer. Funct. Anal. Optim. 30, 613–631 (2009)
 19.
Šamajová, H., Li, T.: Oscillators near Hopf bifurcation. Komunikácie (Žilina) 17, 83–87 (2015)
 20.
Rebenda, J., Šmarda, Z.: A differential transformation approach for solving functional differential equations with multiple delays. Commun. Nonlinear Sci. Numer. Simul. 48, 246–257 (2017)
 21.
Yang, X.J., Tenreiro Machado, J.A., Srivastava, H.M.: A new numerical technique for solving the local fractional diffusion equation: twodimensional extended differential transform approach. Appl. Math. Comput. 274, 143–151 (2016)
 22.
Šamajová, H.: Semianalytical approach to initial problems for systems of nonlinear partial differential equations with constant delay. In: Mikula, K., Sevcovic, D., Urban, J. (eds.) Proceedings of EQUADIFF 2017 Conference, pp. 163–172. Spektrum STU Publishing, Bratislava (2017)
 23.
Rebenda, J., Šmarda, Z., Khan, Y.: A new semianalytical approach for numerical solving of Cauchy problem for differential equations with delay. Filomat 31(15), 4725–4733 (2017)
 24.
Šmarda, Z., Diblík, J., Khan, Y.: Extension of the differential transformation method to nonlinear differential and integrodifferential equations with proportional delays. Adv. Differ. Equ. 2013, 69 (2013)
 25.
Šmarda, Z., Khan, Y.: An efficient computational approach to solving singular initial value problems for Lane–Emden type equations. J. Comput. Appl. Math. 290, 65–73 (2015)
 26.
Rebenda, J.: An application of Bell polynomials in numerical solving of nonlinear differential equations. In: 17th Conference on Applied Mathematics, APLIMAT 2018—Proceedings, pp. 891–900. Spektrum STU, Bratislava (2018)
 27.
Bellen, A., Zennaro, M.: Numerical Methods for Delay Differential Equations. Oxford University Press, Oxford (2003)
 28.
Warne, P.G., Polignone Warne, D.A., Sochacki, J.S., Parker, G.E., Carothers, D.C.: Explicit apriori error bounds and adaptive error control for approximation of nonlinear initial value differential systems. Comput. Math. Appl. 52, 1695–1710 (2006)
 29.
Widatalla, S., Koroma, M.A.: Approximation algorithm for a system of pantograph equations. J. Appl. Math. 9, Article ID 714681 (2012)
 30.
Komashynska, I., AlSmadi, M., AlHabahbeh, A., Ateiwi, A.: Analytical approximate solutions of systems of multipantograph delay differential equations using residual powerseries method. Aust. J. Basic Appl. Sci. 8(10), 664–675 (2014)
Funding
The first author was supported by the Grant CEITEC 2020 (LQ1601) with financial support from the Ministry of Education, Youth and Sports of the Czech Republic under the National Sustainability Programme II. The work of the second author was supported by the Grant FEKTS174225 of Faculty of Electrical Engineering and Communication, Brno University of Technology.
Author information
Affiliations
Contributions
All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Rebenda, J., Šmarda, Z. Numerical algorithm for nonlinear delayed differential systems of nth order. Adv Differ Equ 2019, 26 (2019). https://doi.org/10.1186/s1366201919613
Received:
Accepted:
Published:
MSC
 34K28
 34K07
 34K40
 65L03
Keywords
 Differential transformation
 Method of steps
 Delayed differential system
 Multiple delays