Abstract

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 with ,Now, are :Finally, coefficients for areAt this moment, we substitute and into systems (60) and (61). The next three coefficients at for areand for , we obtainUsing the inverse DT, again we get approximate solution for the IVPs (48), (49), and (50) on the interval :As the calculations are getting more complicated, all the calculations have been done numerically only.

3.4. Numerical Results and Discussion

Table 1 shows comparison of results for Example 1 obtained by DT algorithm with the orders of Taylor polynomials of the approximate solution to results of Matlab function DDENSD in the interval . Since the exact solution is known, absolute errors illustrate precision of each algorithm setting. All numbers are rounded to four decimal places. We see that DDENSD performs satisfactory well and DT for does even better, whereas DT for does not show satisfactory precision.

Table 2 brings the same comparison in the second interval . We can observe a fast growth rate of the function values of the exact solution, which leads to the growth of absolute errors and loss of precision in all settings. It indicates that at the end of the considered interval , the rate of precision would be better seen using relative errors.

Implementation of DT in Matlab in case of Example 2 produces numerical results which are listed in Table 3. The results of DT with order of the Taylor polynomial are compared to values obtained by DT combined with modified Adomian formula in [18] and to values produced by Matlab function DDENSD.

First, we should say that the function DDENSD had difficulty at 0 where the value of the delayed argument was equal to the argument itself. Hence, to make DDENSD work, we replaced by in the second equation of (2). Our hypothesis is that the reason of the DDENSD failure is a combination of two facts: the second equation is neutral with respect to a proportional delay and the interval where the problem is considered contains 0.

Second, we should mention that the numerical results for DDENSD were obtained by looking for approximate solutions on the whole interval . When trying to follow the method of steps, i.e., using DDENSD on and then on , the results on the second interval did not correspond to reality: there was a discontinuity in at .

Furthermore, we recall that the values taken from [18] have been computed using symbolic software Maple and the source code of the computation has been lost.

Now, we can see a very good concordance of all algorithms in numerical values of the second component , while we observe a growing distance between the values of the first component computed by presented DT algorithm and values computed by the other two algorithms. As has exponential characteristics, we interpret the growing distance as growing lack of precision of DT algorithm which is based on approximation by Taylor polynomials. We suppose that dividing the intervals and into smaller subintervals, i.e., refining the mesh grid, and applying the DT algorithm on those smaller intervals consecutively will improve the performance of the presented algorithm.

Although it seems that the algorithm used in [18] shows better performance than the one presented in this paper, we cannot claim it with certainty as the source code got lost and we are not able to reproduce the data. Moreover, the approach used in [18] involves calculations of symbolic derivatives which makes it difficult to implement in numerical software like Matlab.

4. Conclusion

In the paper, we presented an algorithm which makes use of the differential transformation to initial value problems for systems of delayed or neutral differential equations with nonconstant delays. Two examples have been chosen to validate and test the algorithm. Numerical comparison of the presented semianalytical approach to Matlab function DDENSD brought interesting and promising results.

Example 1 showed expected and reliable behaviour of the differential transform in the first step of the method of steps and expected deviation in the numerical results from values of the exact solution in the second step. Furthermore, we could observe a good concordance between the presented algorithm and DDENSD.

After facing difficulties with DDENSD in Example 2, we could confirm a very good concordance of both differential transform and DDENSD in values of the component which has a polynomial character on the considered intervals. On the other hand, we observed a growing discrepancy between the two methods in values of the component which has an exponential character. Our conclusion is that the disagreement is caused by large lengths of the intervals where the approximate solution is computed using the differential transform and that refining the mesh grid is necessary to obtain better performance.

Further investigation will be focused on experimenting with different densities of mesh grids and studying convergence of the algorithm to find the optimal mesh grid. Numerical experiments will be focused on tuning the performance on problems with high complexity whose exact solutions are known and subsequently on applications to nonartificial real-life problems whose exact solutions are unknown.

Data Availability

The numerical data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

The first author was supported by the project CEITEC 2020 (LQ1601) with financial support from the Ministry of Education, Youth and Sports of the Czech Republic under the National Sustainability Programme II, and by the grant FEKT-S-17-4225 of the Faculty of Electrical Engineering and Communication of Brno University of Technology. The Article Processing Charge was funded by Open Access Fund of Brno University of Technology. This support is gratefully acknowledged.

Supplementary Materials

We have included the commented Matlab source code for the first step of the method of steps in Example 1 as supplementary material. (Supplementary Materials)