An algorithm using the differential transformation which is convenient for finding numerical solutions to initial value problems for functional differential equations is proposed in this paper. We focus on retarded equations with delays which in general are functions of the independent variable. The delayed differential equation is turned into an ordinary differential equation using the method of steps. The ordinary differential equation is transformed into a recurrence relation in one variable using the differential transformation. Approximate solution has the form of a Taylor polynomial whose coefficients are determined by solving the recurrence relation. Practical implementation of the presented algorithm is demonstrated in an example of the initial value problem for a differential equation with nonlinear nonconstant delay. A two-dimensional neutral system of higher complexity with constant, nonconstant, and proportional delays has been chosen to show numerical performance of the algorithm. Results are compared against Matlab function DDENSD.

1. Introduction

Functional differential equations (FDEs) are used to model processes and phenomena which depend on past values of the modelled entities. Indicatively, we mention models describing machine tool vibrations [1], predator-prey type models [2], and models used in economics [3]. Further models and details can be found for instance in [4, 5] or [6].

Differential transformation (DT), a semianalytical approach based on Taylor’s theorem, has been proved to be efficient in solving a variety of initial value problems (IVPs), ranging from ordinary to functional, partial, and fractional differential equations [711]. However, there is no publication about systematic application of DT to IVP for differential equations with nonconstant delays which are functions of the independent variable.

In this paper, we present an extension of DT to a class of IVPs for delayed differential equations with analytic righthand side. Albeit the analyticity assumption seems to be quite restrictive, it is reasonable to develop theory for such class of equations [12, 13].

The paper is organised as follows. In Section 2, we define the subject of our study and briefly describe the methods we combine to solve the studied problem, including recalling necessary results of previous studies. Section 3 contains the main results of the paper, including algorithm description, new theorems, examples, and comparison of numerical results. In Section 4, we briefly summarise what has been done in the paper.

2. Methods

2.1. Problem Statement

The problem studied in this paper is to find a solution on a given finite interval to an IVP for the following system of p functional differential equations of n-th order with multiple delays :where , , and are p-dimensional vector functions, are -dimensional vector functions, , , , and are real functions for , where .

We assume that each , where for , , is in general a real function, that is, a time-dependent or time-varying delay. Constant and proportional delays are considered as special cases. In case that some is a proportional delay, we do not require the condition to be valid at 0 if .

Let and ; hence, and . If , we have a retarded system (1); otherwise, if , we call the system neutral. Furthermore, if , initial vector function must be prescribed on the interval .

DT algorithm for the case with all delays being proportional is described in [14]. DT algorithm for the case when all delays are constant is introduced in [15]. In this paper, we develop the algorithm for the case when at least one delay is nonconstant.

To have a complete IVP, we consider system (1) together with initial conditions:and, since , also subject to initial vector function on interval such that

We consider the IVPs (1)–(3) under the following hypotheses:(H1) We assume that all the functions , , are analytic in , the functions , , are analytic in and the functions , , are analytic in an open set containing .(H2) If and in for some and , that is, jth equation is neutral with respect to the proportional delay , we assume that for , . This hypothesis is included since if it is not fulfilled, the existence of unique solution of IVP could be violated.

We note that these assumptions imply that the IVP (1)–(3) has a unique solution in the interval .

2.2. Method of Steps

The basic idea of our approach is to combine DT and the general method of steps. The method of steps enables us to replace the terms including delays with initial vector function and its derivatives. Then, the original IVP for the delayed or neutral system of differential equations is turned into IVP for a system of ordinary differential equations.

For the sake of clarity, we include a simple explanatory example. Suppose that we have a system with three delays, one of each type considered: , , and . We have to distinguish two cases:(a)If , applying the method of steps turns system (1) intowhile(b)If , system (1) is simplified towhereand for . More details on the general method of steps can be found, for instance, in monographs [4] or [6].

Continuation of the method of steps algorithm for equations with constant delays is described in [15]. Briefly summarised, the interval is divided into subintervals , , where and , , are the principal discontinuity points which is the set of points , such that and for , are the minimal roots with odd multiplicity of r equations:

If nonconstant nonproportional delays appear in system (1), the principal set of discontinuity points is defined as follows:

Definition 1. The principal discontinuity points for the solutions of system (1) are given by the set of points , such that and for , are the minimal roots with odd multiplicity of r equations:Similar to the case of constant delays, we break the interval into subintervals , . We start with the mesh grid formed by the principal discontinuity points calculated using Definition 1. To improve convergence or performance of the algorithm, there is a possibility to refine the mesh grid by inserting other points into it. For more details on the principal discontinuity points and mesh grid, we refer to the monograph [16].

2.3. Differential Transformation

Definition 2. Differential transformation of a real function at a point is , where , component of the differential transformation of the function at , , is defined asprovided that the original function is analytic in a neighbourhood of .

Definition 3. Inverse differential transformation of is defined asIn applications, the function is expressed by a finite sumAs we can observe in (10), DT is based on Taylor series; hence, any theorem about convergence of Taylor series may be used. However, we would like to point out the paper [17] where the finest general explicit a priori error estimates are given.
The following formulas are listed, e.g., in [18] and will be used in Section 3.3.

Lemma 1. Assume that and are differential transformations of functions and , respectively:

Remark 1. Similar formulas can be obtained using numerical approach called Functional Analytical Technique based on Operator Theory [19, 20].
The main disadvantage of many papers about DT is that there are almost no examples of equations with nonpolynomial nonlinear terms containing unknown function like, for instance, or . However, DT of components containing nonlinear terms can be obtained in a consistent way using the algorithm described in [21].

Theorem 1. Let and f be real functions analytic near and , respectively, and let h be the composition . Denote , , and as the differential transformations of functions , f, and h at , , and , respectively. Then, the numbers in the sequence satisfy the relations andwhere are the partial ordinary Bell polynomials.

The following Lemma proved in [21] is useful when calculating partial ordinary Bell polynomials.

Lemma 2. The partial ordinary Bell polynomials , , satisfy the recurrence relationwhere and for .

3. Results and Discussion

3.1. Algorithm Description

Recall system (1)with initial conditionsand initial vector function on interval satisfying

Further recall that in Section 2.2, we broke the interval into subintervals , . Define .

Then, we are looking for a solution of the IVP (1)–(3) in the formwhere solution in the jth interval is obtained in the following way. We solve the following equation:where

In case that for , then again

Application of DT at to equation (19) yields a system of recurrence algebraic equations:where the function is the DT of the righthand side of equation (19) and involves application of Theorem 1.

Next, we transform the initial conditions (2). Following Definition 2, we derive

Using (22) with (23) and then inverse transformation rule, we obtain approximate solution to (19) in the form of Taylor series:for all .

To transform (20) correctly, we need the following theorem.

Theorem 2. Let for , where . Let . Denote . Then,where , for , andfor .

Proof. To prove (25) with , we use Theorem 1 with , , and . We immediately getfor . For , Theorem 1 yields . Now, (25) for is a consequence of Lemma 1 and it remains to prove (26). We recall thatAs the assumption was that , we may apply Definition 2 to (28) and obtainSubstituting and gives (26).

3.2. New DT Formulas

In the applications, we also use the following DT formulas.

Theorem 3. Assume that is the differential transformation of the function and :(a)If , then for all t such that ,where and is the Pochhammer symbol.(b)If , then for .

Proof. (a)Recall the Newton’s generalisation of the binomial formula: if x and y are real numbers with , and r is any complex number, one haswhere . Let us rewrite as . Applying (30) yields(b)We start by proving the formulaby induction. For , we have ; hence, (32) is valid. Suppose that (32) holds for k. Then,Thus, formula (32) is valid for all . Now by Definition 2,

3.3. Applications

In this section, we introduce two test problems and show how the practical implementation of the presented algorithm looks like in concrete examples. Comparison of numerical results is given in Section 3.4.

As the first test problem, we choose an IVP for a scalar equation with one nonconstant delay where the exact solution is known to be the exponential function . The purpose of including this example is to compare results obtained by DT against values of the exact solution and also against results obtained by Matlab function DDENSD which has been designed to approximate solutions to IVP for neutral delayed differential equations.

Example 1. Consider the delayed equation:with the initial conditionand with the initial functionFirst we find the differential transform of the initial condition (36) which is . Further denote as the transformation of the exponential function with the center at 0 and as the transformation of the logarithmic function at 1, respectively. Then, Lemma 1 and Theorem 3 yieldFor , equation (35) is transformed intowhereWe haveUsing the inverse transformation, we see that for ,which corresponds to the exact solution to the IVPs (35)–(37).
In the second step of the method of steps, i.e., in the interval , we know that for and equation (35) is transformed intowhereHere, , according to Theorem 3, are coefficients of Taylor series of logarithmic function with the center at :Taking the values calculated in the first step and substituting them into the recurrence formulas (43) and (44), we obtainHence, for , we havewhich again coincides with the exact solution to problems (35)–(37).
In the second application, we have chosen an IVP for a nonlinear system of neutral delayed differential equations taken from the fully open access paper [18]. There are several reasons to test the proposed algorithm on the particular problem. The first is that the problem involves a nonlinear system of neutral equations of high complexity whose exact solution is unknown. Secondly, the proposed algorithm is a complete differential transform version of the algorithm presented in [18] where modified Adomian formula has been used. Furthermore, the calculations done in [18] are shown only for the first step of the method of steps up to the first principal discontinuity point, whereas we continue calculations beyond that point in this paper. Last but not least, we want to verify performance and reproduce values obtained by DT and published in [18]. Rebenda et al. [18] has been submitted 4 years ago for the first time, and since that time, the Maple source code has been lost.

Example 2. Consider a nonlinear system of neutral delayed differential equations:with initial functionsfor , and initial conditionsFor , where is the minimal root of , using the method of steps, we obtainWe need to find the differential transform of the considered problem. We notice that system (2) contains nonlinear term . To get DT of this term, , and we apply Theorem 1. First, applying DT to system (2) at , we get the recurrent system:Denote ; then, , and following Theorem 1, we obtain where and, Theorem 3 being applied, for . Furthermore, the transformed initial conditions areUsing them, we compute the first three coefficients of the nonlinear term :Solving recurrent systems (52) and (53), we getUsing the inverse DT (Definition 3), we get approximate solution for the IVPs (48)–(50) on the interval :which is exactly the same approximate solution which has been obtained in [18].
The second step brings us to solving the given IVP on the interval , where is the minimal root of , . Now taking into account that both proportional delays and and also the time-dependent delay map the interval into the interval , system (2) becomesDenote and , where . By application of Theorem 2 to corresponding terms, system (2) transformed at reads asNow denote ; then, according to Theorem 3, for and Theorem 1 impliesFurther denote and . Then, and, since , Theorem 1 in combination with Lemma 1 yieldsTo get the initial data and for , we have to transformat , . For , we havei.e.,The initial values at t1 will be approximated by taking finite sums in computer evaluations of the infinite sums above. Observe that .
Now let us compute the first few values of . Denote . Then, the first values of areand the coefficients for areLet us turn our attention to the first few values of . Starting