Numerical Algorithm for Nonlinear Delayed Differential Systems of $n$th Order

The purpose of this paper is to propose a semi-analytical technique convenient for numerical approximation of solutions of the initial value problem for $p$-dimensional 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 all three types of delays considered. Accuracy of the results is compared to results obtained by the Laplace decomposition algorithm, the residual power series method and Matlab package DDENSD. Comparison of computing time is done too, showing reliability and efficiency of the proposed technique.


Introduction
Many engineering and chemical processes as well as economic and biological systems are modeled by a set of nonlinear delay differential equations (FDEs). As examples, we mention models describing machine tool vibrations [1], behaviour of the central nervous system in a learning process, species populations struggling for a common food, dynamics of an autogenerator with delay and second-order filter, evolution of population of one or more species etc. Further models and details can be found for instance in [2] or [3].
To study the local dynamics of FDEs, the nonlinearities are expanded in a Taylor series about a fixed point or a periodic solution. The latter case gives rise to linear and nonlinear FDEs, depending on the degree of the series, with time-periodic coefficients. Recent example of Taylor series approach applied to study the local dynamics can be found in [4].
In the last two decades, various methods such as homotopy analysis method (HAM), homotopy perturbation method (HPM), variational iteration method (VIM), Adomian decomposition method (ADM), and also variety of methods derived from Taylor series method such as Taylor polynomial method, Taylor collocation method and differential transformation method (DTM) have been considered to approximate solutions of certain classes of FDEs in a series form. However, in several papers, for example Evans and Raslan [5], initial problems are not properly defined. The authors used only initial conditions in certain points, not the initial function on the whole interval, thus the way to obtain solutions of illustrative examples is not correct. Moreover, formulas used in the calculations are very complicated. Similar situation occured in paper Blanco-Cocom, Estrella and Avila-Vales [6], where the algorithms are not able to approximate a solution starting with another initial function. Both mentioned papers are using ADM. Comparison of ADM and DTM, which we decided to employ, is done in [7].
The concept of differential transformation in the form we utilize was proposed by Zhou [8] in 1986. He applied the method to solve linear and nonlinear initial value problems in electric circuit analysis. This method constructs a semi-analytical numerical technique that uses Taylor method for solving of differential equations in the form of a polynomial. It is different from high-order Taylor series method which requires symbolic computation of necessary derivatives of the data functions. This is not the only way how to employ Taylor approach to solving FDEs. One possibility is to start with characterictic quasipolynomials and proceed to polynomial quasisolutions as in [9] and [10]. It is also possible to derive the formulas using more general approach involving Hilbert spaces, see [11] or [12].
The crucial idea of our Taylor series concept is to combine general method of steps suitable for Cauchy problems for FDEs with a Taylor series technique known as DTM. For more details on method of steps see e.g. monographs Kolmanovskii and Myshkis [2], Hale and Verduyn Lunel [3] or Bellen and Zennaro [13]. Our approach enables us to replace the terms involving delay with initial function and its derivatives. Consequently, the original Cauchy problem for delayed or neutral differential equation is reduced to Cauchy problem for ordinary differential equation. Also the ambiguities mentioned above are removed. Moreover, while ADM, HAM, HPM and VIM require initial approximation guess and symbolic computation of necessary derivatives and, in general, n-dimensional integrals in iterative schemes, presented algorithm is different: Cauchy problem for FDE is reduced to a system of recurrent algebraic relations.

Differential Transformation
Differential transformation of the k-th derivative of function u(t) is defined as follows: where u(t) is the original function and U(k) is the transformed function. Inverse differential transformation of U(k) is defined as follows: hence In fact, formula (6) indicates that the concept of differential transformation is derived from Taylor series expansion. However, DTM is not evaluating the symbolic derivatives: instead of this, relative derivatives are calculated by an iterative way which is described by the transformed form of the original equation.

Lemma 1. Assume that F(k), H(k) and U(k) are differential transformations of functions f (t), h(t) and u(t), respectively. Then:
If f (t) = t n , then F(k) = δ(k − n), where δis the Kronecker delta. If Similarly transformation formulas for proportionally delayed arguments can be proved, as in [14]: Assume that W(k), U(k), U i (k) are the differential transformations of the functions w(t), u(t), u i (t) and q, q i ∈ (0, 1), i = 1, 2 then: The main disadvantage of standard form of DTM is that it cannot be used directly for equations with nonlinear terms containing unknown function u(t). For example, for f (u) = √ 1 + u 2 or f (u) = e sin u etc. However, differential transformation of components containing nonlinear terms can be calculated using the socalled Adomian polynomials A n in which each solution component u i is replaced by the corresponding differential transformations component U(i), i = 0, 1, 2, . . . (see [15]). If F(k) is the differential transformation of a nonlinear term f (u), then For illustration, we give several components explicitly: . . .
Next we transform initial conditions (2). Following definition of differential transformation, we derive Applying DTM to system (8) we get system of recurrent algebraic equations Using transformed initial conditions and then inverse transformation rule, we obtain approximate solution of system (1) in the form of infinite Taylor series If t * = 0, which means that all delayed arguments are proportional, the approximate solution u(t) is valid on the intersection of its convergence interval with [0, T * ]. Otherwise, if t * < 0, we denote t α i = inf{t : α i (t) > 0} and t α = min  (2) and (3) eventually, is locally approximable by using combination of method of steps and differential transformation method or differential transformation only.
Proof. If all the assumptions are fulfilled, then either we are dealing with the system of equations with proportional delays only, or we can pass to such system or system of ODEs using the method of steps. Validity of this approach can be verified for instance in [2], [3] or [13].
The obtained system has a unique solution subject to initial conditions (2), either on [0, T * ] or on [0, t α ] ∩ [0, T * ]. Then, since the righthand side of the system is analytic, differential transformation leads to solving a system of difference equations instead of solving the system of differential equations.
Using the transformed initial conditions, a finite part of the sequence which is a solution of difference equation is calculated. From the definition of differential transform we see that the finite sequence we calculated represents coefficients of Taylor polynomial of the solution of Cauchy problem. According to the Theorem 2.2.1 in [13], this polynomial is a local approximation of the unique solution which is valid on some interval [0, δ], 0 < δ ≤ T * α where T * α = min{t α , T * }.
where N is the order of Taylor polynomial approximation and K is a bound for (N + 1)-th derivative of the solution on [0, δ].

Proof.
Since the functions f j are analytical in [0, T * ] × R np × R ωp , without the loss of generality we can suppose that using Theorem 1 we are able to approximate the unique solution by Taylor polynomial The remainder of Taylor series is given by for some nonnegative constant K, then the maximum error for u N (t) in this interval can be estimated from this remainder term as e N = K (N + 1)! .
Example 3. Consider the nonlinear system of neutral delayed differential equations with initial functions for t ∈ [−2, 0] and initial conditions Using the method of steps we get System (22) can not be solved using classic DTM with respect to nonlinear term f (u) = 3 u 2 1 , therefore we will apply modified Adomian formula for differential transformation components to nonlinear term f (u) . At first, applying DTM to system (22) we get recurrent system (k + 1)(k + 2)(k + 3)U 1 (k + 3) = e −2 k l=0 1 l!
• There is no need for calculating multiple integrals or derivatives and less computational work is demanded compared to other popular methods (Adomian decomposition method, variational iteration method, homotopy perturbation method, homotopy analysis method).
• Using presented approach, we are able not only to obtain approximate solution, but even there is a possibility to identify unique solution of Cauchy problem in closed form.
• A specific advantage of this technique over any purely numerical method is that it offers a smooth, functional form of the solution over a time step.
• Another advantage is that using our approach we avoided ambiguities, incorrect formulations and ill-posed problems that occured in recent papers.
• Finally, a subject of further investigation is to develop the presented technique for system (1) with state dependent delays and for partial differential equations.