# Differential Transform Algorithm for Functional Differential Equations with Time-Dependent Delays

Authors: Josef Rebenda; Zuzana Pátíková
Journal: Complexity
 (2020)

---

## 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.


---

## Body

## 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 [7–11]. 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 Section2, 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 intervalt0,T⊂0,∞ to an IVP for the following system of p functional differential equations of n-th order with multiple delays α1t,…,αrt:(1)unt=ft,ut,u′t,…,un−1t,u1α1t,…,urαrt,where unt=u1nt,…,upntT,ukt=u1kt,…,upkt, k=0,1,…,n−1, and f=f1,…,fpT are p-dimensional vector functions, uiαit=uαit,u′αit,…,umiαit are mi⋅p-dimensional vector functions, mi≤n, i=1,2,…,r, r∈ℕ, and fj:t0,T are real functions for j=1,2,…,p, where ω=∑i=1rmi.We assume that eachαit=t−τit, where τit≥τi0>0 for t∈t0,T, i=1,2,…,r, 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 αi is a proportional delay, we do not require the condition τit≥τi0>0 to be valid at 0 if t0=0.Lett∗=min1≤i≤rinft∈t0,Tαit and m=maxm1,m2,…,mr; hence, t∗≤t0 and m≤n. If m<n, we have a retarded system (1); otherwise, if m=n, we call the system neutral. Furthermore, if t∗<0, initial vector function Φt=ϕ1t,…,ϕptT must be prescribed on the interval t∗,t0.DT algorithm for the caset∗=t0=0 with all delays being proportional is described in [14]. DT algorithm for the case t∗<t0 when all delays are constant is introduced in [15]. In this paper, we develop the algorithm for the case t∗<t0 when at least one delay is nonconstant.To have a complete IVP, we consider system (1) together with initial conditions:(2)ut0=v0,u′t0=v1,…,un−1t0=vn−1,and, since t∗<t0, also subject to initial vector function Φt on interval t∗,t0 such that(3)Φt0=ut0,…,Φn−1t0=un−1t0.We consider the IVPs (1)–(3) under the following hypotheses:(H1) We assume that all the functionsϕjt, j=1,…,p, are analytic in t∗,t0, the functions αit, i=1,…,r, are analytic in t0,T and the functions fj, j=1,…,p, are analytic in an open set containing t0,T×ut0,uT×…×urαrt0,urαrT.(H2) Ifαit=qit and mi=n in fj for some i∈1,…,r and j∈1,…,p, that is, jth equation is neutral with respect to the proportional delay αi, we assume that ulnαit≡0 for l∈1,…,p, l≠j. 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 t0,T.


### 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Φt 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:α1t=t−τ1t, α2t=t−τ2, and α3t=q3t. We have to distinguish two cases:(a)
Ift0=0, applying the method of steps turns system (1) into(4)unt=ft,ut,…,un−1t,Φ1t−τ1t,Φ2t−τ2,u3q3t,while(b)
Ift0>0, system (1) is simplified to(5)unt=ft,ut,…,un−1t,Φ1t−τ1t,Φ2t−τ2,Φ3q3t,where(6)Φ1t−τ1t=Φt−τ1t,Φ′t−τ1t,…,Φm1t−τ1t,Φ2t−τ2=Φt−τ2,Φ′t−τ2,…,Φm2t−τ2,u3q3t=uq3t,u′q3t,…,um3q3t,Φ3q3t=Φq3t,Φ′q3t,…,Φm3q3t,and ml≤n for l=1,2,3,4. 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τ1,…,τr is described in [15]. Briefly summarised, the interval t0,T is divided into subintervals Il=tl−1,tl, l=1,…,K, where tK=T and tl, l=1,…,K−1, are the principal discontinuity points which is the set of points tρ,σ, such that t0,1=t0 and for ρ,σ≥1, tρ,σ are the minimal roots with odd multiplicity of r equations:(7)tρ,σ−1r+μ−τμ=tρ−1,σ,μ=1,…,r.If nonconstant nonproportional delaysαi 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 tρ,σ, such that t0,1=t0 and for ρ,σ≥1, tρ,σ are the minimal roots with odd multiplicity of r equations:(8)αμtρ,σ−1r+μ=tρ−1,σ,μ=1,…,r.
Similar to the case of constant delays, we break the intervalt0,T into subintervals Il=tl−1,tl, l=1,…,K. We start with the mesh grid t0,…,tK 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 functionut at a point t0∈ℝ is Dutt0=Ukt0k=0∞, where Ukt0, k−th component of the differential transformation of the function ut at t0, k∈ℕ0, is defined as(9)Ukt0=1k!dkutdtkt=t0,provided that the original function ut is analytic in a neighbourhood of t0.Definition 3.
Inverse differential transformation ofUkt0k=0∞ is defined as(10)ut=D−1Ukt0k=0∞t0=∑k=0∞Ukt0t−t0k.
In applications, the functionut is expressed by a finite sum(11)ut=∑k=0NUkt0t−t0k.
As 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 thatFkt0 and Ukt0 are differential transformations of functions ft and ut, respectively:(12)Ifft=dnutdtn,thenFkt0=k+n!k!Uk+nt0.Ifft=tn,thenFk0=δk−n,whereδk−n=δknistheKroneckerdelta.Ifft=eλt,thenFk0=λkk!.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 functionut like, for instance, fu=1+u35 or fu=esinu. However, DT of components containing nonlinear terms can be obtained in a consistent way using the algorithm described in [21].Theorem 1.
Letg and f be real functions analytic near t0 and gt0, respectively, and let h be the composition ht=f∘gt=fgt. Denote Dgtt0=Gkk=0∞, Dftgt0=Fkk=0∞, and Df∘gtt0=Hkk=0∞ as the differential transformations of functions g, f, and h at t0, gt0, and t0, respectively. Then, the numbers Hk in the sequence Hkk=0∞ satisfy the relations H0=F0 and(13)Hk=∑l=1kFl⋅B^k,lG1,…,Gk−l+1,fork≥1,where B^k,lx^1,…,x^k−l+1 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 polynomialsB^k,lx^1,…,x^k−l+1, l=1,2,…,k≥l, satisfy the recurrence relation(14)B^k,lx^1,…,x^k−l+1=∑i=1k−l+1i⋅lkx^iB^k−i,l−1x^1,…,x^k−i−l+2,where B^0,0=1 and B^k,0=0 for k≥1.


## 2.1. Problem Statement

The problem studied in this paper is to find a solution on a given finite intervalt0,T⊂0,∞ to an IVP for the following system of p functional differential equations of n-th order with multiple delays α1t,…,αrt:(1)unt=ft,ut,u′t,…,un−1t,u1α1t,…,urαrt,where unt=u1nt,…,upntT,ukt=u1kt,…,upkt, k=0,1,…,n−1, and f=f1,…,fpT are p-dimensional vector functions, uiαit=uαit,u′αit,…,umiαit are mi⋅p-dimensional vector functions, mi≤n, i=1,2,…,r, r∈ℕ, and fj:t0,T are real functions for j=1,2,…,p, where ω=∑i=1rmi.We assume that eachαit=t−τit, where τit≥τi0>0 for t∈t0,T, i=1,2,…,r, 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 αi is a proportional delay, we do not require the condition τit≥τi0>0 to be valid at 0 if t0=0.Lett∗=min1≤i≤rinft∈t0,Tαit and m=maxm1,m2,…,mr; hence, t∗≤t0 and m≤n. If m<n, we have a retarded system (1); otherwise, if m=n, we call the system neutral. Furthermore, if t∗<0, initial vector function Φt=ϕ1t,…,ϕptT must be prescribed on the interval t∗,t0.DT algorithm for the caset∗=t0=0 with all delays being proportional is described in [14]. DT algorithm for the case t∗<t0 when all delays are constant is introduced in [15]. In this paper, we develop the algorithm for the case t∗<t0 when at least one delay is nonconstant.To have a complete IVP, we consider system (1) together with initial conditions:(2)ut0=v0,u′t0=v1,…,un−1t0=vn−1,and, since t∗<t0, also subject to initial vector function Φt on interval t∗,t0 such that(3)Φt0=ut0,…,Φn−1t0=un−1t0.We consider the IVPs (1)–(3) under the following hypotheses:(H1) We assume that all the functionsϕjt, j=1,…,p, are analytic in t∗,t0, the functions αit, i=1,…,r, are analytic in t0,T and the functions fj, j=1,…,p, are analytic in an open set containing t0,T×ut0,uT×…×urαrt0,urαrT.(H2) Ifαit=qit and mi=n in fj for some i∈1,…,r and j∈1,…,p, that is, jth equation is neutral with respect to the proportional delay αi, we assume that ulnαit≡0 for l∈1,…,p, l≠j. 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 t0,T.


## 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Φt 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:α1t=t−τ1t, α2t=t−τ2, and α3t=q3t. We have to distinguish two cases:(a)
Ift0=0, applying the method of steps turns system (1) into(4)unt=ft,ut,…,un−1t,Φ1t−τ1t,Φ2t−τ2,u3q3t,while(b)
Ift0>0, system (1) is simplified to(5)unt=ft,ut,…,un−1t,Φ1t−τ1t,Φ2t−τ2,Φ3q3t,where(6)Φ1t−τ1t=Φt−τ1t,Φ′t−τ1t,…,Φm1t−τ1t,Φ2t−τ2=Φt−τ2,Φ′t−τ2,…,Φm2t−τ2,u3q3t=uq3t,u′q3t,…,um3q3t,Φ3q3t=Φq3t,Φ′q3t,…,Φm3q3t,and ml≤n for l=1,2,3,4. 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τ1,…,τr is described in [15]. Briefly summarised, the interval t0,T is divided into subintervals Il=tl−1,tl, l=1,…,K, where tK=T and tl, l=1,…,K−1, are the principal discontinuity points which is the set of points tρ,σ, such that t0,1=t0 and for ρ,σ≥1, tρ,σ are the minimal roots with odd multiplicity of r equations:(7)tρ,σ−1r+μ−τμ=tρ−1,σ,μ=1,…,r.If nonconstant nonproportional delaysαi 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 tρ,σ, such that t0,1=t0 and for ρ,σ≥1, tρ,σ are the minimal roots with odd multiplicity of r equations:(8)αμtρ,σ−1r+μ=tρ−1,σ,μ=1,…,r.
Similar to the case of constant delays, we break the intervalt0,T into subintervals Il=tl−1,tl, l=1,…,K. We start with the mesh grid t0,…,tK 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 functionut at a point t0∈ℝ is Dutt0=Ukt0k=0∞, where Ukt0, k−th component of the differential transformation of the function ut at t0, k∈ℕ0, is defined as(9)Ukt0=1k!dkutdtkt=t0,provided that the original function ut is analytic in a neighbourhood of t0.Definition 3.
Inverse differential transformation ofUkt0k=0∞ is defined as(10)ut=D−1Ukt0k=0∞t0=∑k=0∞Ukt0t−t0k.
In applications, the functionut is expressed by a finite sum(11)ut=∑k=0NUkt0t−t0k.
As 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 thatFkt0 and Ukt0 are differential transformations of functions ft and ut, respectively:(12)Ifft=dnutdtn,thenFkt0=k+n!k!Uk+nt0.Ifft=tn,thenFk0=δk−n,whereδk−n=δknistheKroneckerdelta.Ifft=eλt,thenFk0=λkk!.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 functionut like, for instance, fu=1+u35 or fu=esinu. However, DT of components containing nonlinear terms can be obtained in a consistent way using the algorithm described in [21].Theorem 1.
Letg and f be real functions analytic near t0 and gt0, respectively, and let h be the composition ht=f∘gt=fgt. Denote Dgtt0=Gkk=0∞, Dftgt0=Fkk=0∞, and Df∘gtt0=Hkk=0∞ as the differential transformations of functions g, f, and h at t0, gt0, and t0, respectively. Then, the numbers Hk in the sequence Hkk=0∞ satisfy the relations H0=F0 and(13)Hk=∑l=1kFl⋅B^k,lG1,…,Gk−l+1,fork≥1,where B^k,lx^1,…,x^k−l+1 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 polynomialsB^k,lx^1,…,x^k−l+1, l=1,2,…,k≥l, satisfy the recurrence relation(14)B^k,lx^1,…,x^k−l+1=∑i=1k−l+1i⋅lkx^iB^k−i,l−1x^1,…,x^k−i−l+2,where B^0,0=1 and B^k,0=0 for k≥1.


## 3. Results and Discussion

### 3.1. Algorithm Description

Recall system (1)(15)unt=ft,ut,u′t,…,un−1t,u1α1t,…,urαrt,with initial conditions(16)ut0=v0,u′t0=v1,…,un−1t0=vn−1,and initial vector function Φt on interval t∗,t0 satisfying(17)Φt0=ut0,…,Φn−1t0=un−1t0.Further recall that in Section2.2, we broke the interval t0,T into subintervals Il=tl−1,tl, l=1,…,K. Define I0=t∗,t0.Then, we are looking for a solutionut of the IVP (1)–(3) in the form(18)ut=uI1t,t∈I1,uI2t,t∈I2,⋮uIKt,t∈IK,where solution uIj in the jth interval Ij is obtained in the following way. We solve the following equation:(19)uIjnt=ft,uIjt,uIj′t,…,uIjn−1t,uj,1α1t,…,uj,rαrt,where(20)uj,iαit=uIlαit,uIl′αit,…,uIlmiαit,if
α
i
t
∈
I
l
,
for
t
∈
I
j
,
l
∈
1
,
…
,
j
,
j
∈
1
,
…
.
K
.In case thatαit∈I0=t∗,t0 for t∈Ij, then again(21)uj,iαit=ϕαit,ϕ′αit,…,ϕmiαit.Application of DT attj−1 to equation (19) yields a system of recurrence algebraic equations:(22)UIjk+ntj−1=Fk,UIjk,UIjk+1,…,UIjk+n−1,where the function F 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(23)UIjktj−1=1k!uIjktj−1,fork=0,1,…,n−1,j∈1,…,K.Using (22) with (23) and then inverse transformation rule, we obtain approximate solution to (19) in the form of Taylor series:(24)uIjt=∑k=0∞UIjktj−1t−tj−1k,t∈Ij,for all j∈1,…,K.To transform (20) correctly, we need the following theorem.Theorem 2.
Letαit∈Il for t∈Ij, where l∈1,…,j−1. Let p∈ℕ. Denote Dαittj−1=Aiktj−1k=0∞. Then,(25)DuIlpαittj−1=∑y=0ky+p!y!UIly+pαitj−1⋅B^k,yAi1tj−1,…,Aik−y+1tj−1k=0∞,where B^0,0=1, B^k,0=0 for k≥1, and(26)UIlyαitj−1=∑x=0∞x+yxαitj−1−tl−1xUIlx+ytl−1,for y≥0.Proof.
To prove (25) with p=0, we use Theorem 1 with ft=uIlt, gt=αit, and ht=f∘gt. We immediately get(27)Hktj−1=∑y=1kUIlyαitj−1⋅B^k,yAi1tj−1,…,Aik−y+1tj−1,for k≥1. For k=0, Theorem 1 yields H0tj−1=UIl0αitj−1=UIl0αitj−1⋅B^0,0Ai1tj−1. Now, (25) for p>0 is a consequence of Lemma 1 and it remains to prove (26). We recall that(28)uIlt=∑k=0∞UIlktl−1t−tl−1k,t∈Il.
As the assumption was thatαitj−1∈Il, we may apply Definition 2 to (28) and obtain(29)UIlyαitj−1=1y!dyuIltdtyt=t01y!∑z=y∞z!z−y!UIlztl−1t0−tl−1z−y.
Substitutingt0=αitj−1 and z=x+y gives (26).


### 3.2. New DT Formulas

In the applications, we also use the following DT formulas.Theorem 3.
Assume thatFkt0 is the differential transformation of the function ft and r∈R:(a)
Ifft=tr, then Fkt0=rkt0r−k for all t such that t−t0<t0,whererk=rr−1…r−k+1/k!=rk/k! and rk is the Pochhammer symbol.(b)
Ifft=lnt, then Fkt0=−1k−1/k⋅t0k for k≥1.Proof.
(a)
Recall the Newton’s generalisation of the binomial formula: ifx and y are real numbers with x>y, and r is any complex number, one has(30)x+yr=∑k=0∞rkxr−kyk,where rk=rr−1…r−k+1/k!. Let us rewrite tr as tr=t−t0+t0r=t0+t−t0r. Applying (30) yields(31)tr=∑k=0∞rkt0r−kt−t0k.(b)
We start by proving the formula(32)lntk=−1k−1k−1!tk,by induction. For k=1, we have lnt′=1/t; hence, (32) is valid. Suppose that (32) holds for k. Then,(33)lntk+1=lntk′=−1k−1k−1!tk′=−1k−1k−1!t−k′=−1k−1k−1!−kt−k−1=−1kk!tk+1.Thus, formula (32) is valid for all k∈ℕ. Now by Definition 2,(34)Fkt0=1k!dklntdtkt=t0=1k!−1k−1k−1!tkt=t0=−1k−1k⋅t0k.


### 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 Section3.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 functionet. 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:(35)u′t=ut−t+ulnt,with the initial condition(36)u1=e,and with the initial function(37)ϕt=et,t∈0,1.
First we find the differential transform of the initial condition (36) which is U01=e. Further denote Det0=Ek0k=0∞ as the transformation of the exponential function with the center at 0 and Dlnt1=Fk1k=0∞ as the transformation of the logarithmic function at 1, respectively. Then, Lemma 1 and Theorem 3 yield(38)Ek0=1k!,Fk1=−1k−1k,fork≥1.
Fort∈1,e, equation (35) is transformed into(39)k+1Uk+11=Uk1−δk−δk−1+Hk1,where(40)Hk1=∑l=1kEl0B^k,lF11,…,Fk−l+11,fork≥1,H01=E00=1.
We have(41)U11=e−1−0+1=e,H11=E10⋅F11=1⋅1=1,U21=12U11−0−1+H11=12e−1+1=e⋅12,H21=E10⋅F21+E20⋅F112=1⋅−12+12⋅1=0,U31=13U21+H21=13⋅e⋅12=e⋅13!,⋮
Using the inverse transformation, we see that fort∈1,e,(42)ut=e1+t−1+t−122+t−133!…=e⋅et−1=et,which corresponds to the exact solution to the IVPs (35)–(37).
In the second step of the method of steps, i.e., in the intervalt∈e,ee, we know that ut=et for t∈1,e and equation (35) is transformed into(43)k+1Uk+1e=Uke−eδk−δk−1+Hke,where(44)H0e=U01,Hke=∑l=1kUl1B^k,lF1e,…,Fk−l+1e,fork≥1.
Here,Fke, according to Theorem 3, are coefficients of Taylor series of logarithmic function with the center at e:(45)F0e=1,Fke=−1k−1k⋅ek,fork≥1.
Taking the values calculated in the first step and substituting them into the recurrence formulas (43) and (44), we obtain(46)U0e=ue=ee,U1e=U0e−eδ0−δ−1+H0e=ee−e−0+e=ee,H1e=U11⋅F1e=e⋅1e=1,U2e=12U1e−0−1+H1e=12ee−1+1=12ee,H2e=U11⋅B^2,1F1e,F2e+U21⋅B^2,2F1e=e⋅F2e+e2⋅F1e2=e⋅12−1e2+e2⋅1e2=0,U3e=13U2e+H2e=13!ee,H3e=U11⋅B^3,1F1e,F2e,F3e+U21⋅B^3,2F1e,F2e+U31⋅B^3,3F1e=e⋅F3e+e2⋅2⋅F1e⋅F2e+e3!F1e3=e⋅2e313!+e⋅1e⋅−1e2⋅12+e6⋅1e3=1e213−12+16=0,U4e=14U3e+H3e=14!ee.Hence, for t∈e,ee, we have(47)ut=ee+eet−e+12eet−e2+13!eet−e3+14!eet−e4+…=ee1+t−e+t−e22+t−e33!+t−e44!+…=ee⋅et−e=et,which 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:(48)u1‴=u1‴t−2u1t3+u1t23+u2′t−12e−t,u2‴=12u2‴t2+u2′t−1u1t3,with initial functions(49)ϕ1t=et,ϕ2t=t2,for t∈−2,0, and initial conditions(50)u10=1,u2′0=1,u1″0=1,u20=0,u2′0=0,u2″0=2.
Fort∈0,t1, where t1≈0,351734 is the minimal root of t−1/2e−t=0, using the method of steps, we obtain(51)u1‴=et−2u1t3+u1t23+2t−e−t,u2‴=12u2‴t2+2t−1u1t3.
We need to find the differential transform of the considered problem. We notice that system (2) contains nonlinear term ht=u1t23. To get DT of this term, Dht0=H1k0k=0∞, and we apply Theorem 1. First, applying DT to system (2) at t0=0, we get the recurrent system:(52)k+1k+2k+3U1k+30=e−2∑l=0k1l!13k−lU1k−l0+H1k0+2δk−−1kk!,(53)k+1k+2k+31−12k+1U2k+30=213k−1U1k−10−23kU1k0.
Denotegt=t2/3; then, ht=g∘u1t, and following Theorem 1, we obtain(54)H100=G101,H1k0=∑l=1kG1l1B^k,lU110,…,U1k−l+10,fork≥1, where Dgt1=G1k1k=0∞ and, Theorem 3 being applied, G1k1=2/3k for k≥0. Furthermore, the transformed initial conditions are(55)U100=1,U110=1,U120=12,U200=0,U210=0,U220=1.
Using them, we compute the first three coefficients of the nonlinear termht:(56)H100=G101=1,H110=G111⋅B1,1U110=23⋅1=23,H120=G111⋅B2,1U110,U120+G121⋅B2,2U110=23⋅12−19⋅1=29.
Solving recurrent systems (52) and (53), we get(57)k=0:U130=16e−2U100+H100+1=2+e−26,U230=13−2U100=−23.k=1:U140=e−22413U110+U100+124H110+1=4e−2+572,U240=1182U100−23U110=227.k=2:U150=e−26019U120+13U110+12U100+160H120−12=16e−2−51080,U250=210523U110−29U120=2189.
Using the inverse DT (Definition3), we get approximate solution for the IVPs (48)–(50) on the interval 0,t1:(58)u1,I1t=1+t+12t2+2+e−26t3+4e−2+572t4+16e−2−51080t5+…,u2,I1t=t2−23t3+227t4+2189t5+…,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 intervalt1,t2, where t2 is the minimal root of t−1/2e−t=t1, t2≈0,620556. Now taking into account that both proportional delays q1t=1/3t and q2t=1/2t and also the time-dependent delay t−τ1t=t−1/2e−t map the interval t1,t2 into the interval 0,t1, system (2) becomes(59)u1‴=et−2u1,I1t3+u1t23+u2,I1′t−12e−t,u2‴=12u2,I1‴t2+2t−1u1,I1t3.
DenoteDhtt1=H2kt1k=0∞ and Dftt1=F2kt1k=0∞, where ft=u2′t−1/2e−t. By application of Theorem 2 to corresponding terms, system (2) transformed at t0=t1 reads as(60)k+1k+2k+3U1k+3t1=e−2∑l=0ket1l!13k−l∑x=0∞x+k−lxU1x+k−l0t13x+H2kt1+F2kt1,(61)k+1k+2k+3U2k+3t1=12k+1k+3k+2k+1∑x=0∞x+k+3xU2x+k+30t12x+2∑l=0kt1−1δl+δl−113k−l∑x=0∞x+k−lxU1x+k−l0t13x.
Now denoteDgtu1t1=G2ku1t1k=0∞; then, according to Theorem 3, G2ku1t1=2/3ku1t12/3−k for k≥0 and Theorem 1 implies(62)H20t1=G20u1t1=u1t123,H2kt1=∑l=1kG2lu1t1⋅B^k,lU11t1,…,U1k−l+1t1,fork≥1.
Further denoteet=t−1/2e−t and Dett1=E2kt1k=0∞. Then, ft=u2′∘et and, since et1=0, Theorem 1 in combination with Lemma 1 yields(63)E2kt1=t1⋅δk+δk−1−12⋅e−t1−1kk!,k≥0,F20t1=U210=0,F2kt1=∑l=1kl+1U2l+10⋅B^k,lE21t1,…,E2k−l+1t1,fork≥1.
To get the initial dataU1kt1 and U2kt1 for k=0,1,2, we have to transform(64)uit=Ui00+Ui10t+Ui20t2+Ui30t3+…,at t1, i=1,2. For k=0,1,2, we have(65)Uikt1=Ui00+Ui101kt11−k+Ui202kt12−k+Ui303kt13−k+…,i.e.,(66)Ui0t1=∑k=0∞Uik0t1k,Ui1t1=∑k=0∞k+1Uik+10t1k,Ui2t1=∑k=0∞k+22Uik+20t1k.
The initial values att1 will be approximated by taking finite sums in computer evaluations of the infinite sums above. Observe that u1t1=U10t1.
Now let us compute the first few values ofH2kt1. Denote B^k,l=B^k,lU11t1,…,U1k−l+1t1. Then, the first values of B^k,l are(67)B^0,0=1,B^1,1=U11t1⋅B^0,0=U11t1,B^2,1=12U11t1⋅B^1,0+U12t1⋅B^0,0=U12t1,B^2,2=U11t1⋅B^1,1=U11t12,and the coefficients H2 for k=0,1,2 are(68)H20t1=u1t1,H21t1=23u1t1−1/3⋅B^1,1=23u1t1−1/3⋅U11t1,H22t1=23u1t1−1/3⋅B^2,1−19u1t1−4/3⋅B^2,2=23u1t1−1/3⋅U12t1−19u1t1−4/3⋅U11t12.
Let us turn our attention to the first few values ofF2kt1. Starting with E2kt1,(69)E20t1=t1−12e−t1=0,E21t1=1+12e−t1=1+t1,E22t1=−12⋅e−t12=−t12.
Now,B^k,l are B^k,lE21t1,…,E2k−l+1t1:(70)B^0,0=1,B^1,1=E21⋅B^0,0=1+t1,B^2,1=12E21⋅B^1,0+E22⋅B^0,0=−t12,B^2,2=E21⋅B^1,1=1+t12.
Finally, coefficientsF2 for k=0,1,2 are(71)F20t1=U210=0,F21t1=2U220⋅B^1,1=21+t1,F22t1=2U220⋅B^2,1+3U230⋅B^2,2=−t1−21+t12.
At this moment, we substituteH2 and F2 into systems (60) and (61). The next three coefficients at t1 for U1 are(72)k=0:U13t1=16e−2et1⋅u1t13+u1t1,k=1:U14t1=124e−2et1⋅13⋅∑x=0∞x+1U1x+10t13x+et1⋅u1t13+23⋅u1t1−1/3⋅U11t1+21+t1,k=2:U15t1=160e−2et1132⋅∑x=0∞x+2xU1x+20t13x+et1⋅13⋅∑x=0∞x+1U1x+10t13x+et12⋅u1t13+16023u1t1−1/3⋅U12t1−19u1t1−4/3⋅U11t12−t1−21+t12,and for U2, we obtain(73)k=0:U23t1=12∑x=0∞x+3xU2x+30t12x+13t1−1⋅u1t13,k=1:U24t1=14∑x=0∞x+4xU2x+40t12x+112t1−1⋅13⋅∑x=0∞x+1U1x+10t13x+u1t13,k=2:U15t1=18∑x=0∞x+5xU2x+50t12x+130t1−1⋅19⋅∑x=0∞x+2xU1x+20t13x+13∑x=0∞x+1U1x+10t13x⋅
Using the inverse DT, again we get approximate solution for the IVPs (48), (49), and (50) on the intervalt1,t2:(74)u1,I2t=U10t1+U11t1t−t1+U12t1t−t12+U13t1t−t13+…,u2,I2t=U20t1+U21t1t−t1+U22t1t−t12+U23t1t−t13+….
As the calculations are getting more complicated, all the calculations have been done numerically only.


### 3.4. Numerical Results and Discussion

Table1 shows comparison of results for Example 1 obtained by DT algorithm with the orders of Taylor polynomials of the approximate solution N=5,10,25 to results of Matlab function DDENSD in the interval 1,e. 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 N=10,25 does even better, whereas DT for N=5 does not show satisfactory precision.Table 1
Example1, error analysis of u in 1,e.
Exact solution
DT 5
DT 10
DT 25
Matlab DDENSD
t
e
t
u
−
e
t
u
−
e
t
u
−
e
t
u
−
e
t
1
2.7183
0
0
0
0
1.4296
4.1769
0.0000
6.4890E − 12
8.8818E − 16
1.6976E − 5
1.8591
6.4182
0.0017
1.0381E − 8
8.8818E − 16
3.6095E − 5
2.2887
9.8622
0.0211
1.2409E − 6
3.5527E − 15
9.3438E − 5
e
15.1543
0.1273
3.0578E − 5
1.2861E − 12
1.4012E − 4Table2 brings the same comparison in the second interval e,ee. 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 e,ee, the rate of precision would be better seen using relative errors.Table 2
Example1, error analysis of u in e,ee.
Exact solution
DT 5
DT 10
DT 25
Matlab DDENSD
t
e
t
u
−
e
t
u
−
e
t
u
−
e
t
u
−
e
t
e
15.1543
3.5527E − 15
3.5527E − 15
3.5527E − 15
0
5.8273
339.4331
32.3491
1.333
2.2737E − 13
0.0743
8.9363
7.6028E + 3
4.4755E + 3
397.5509
2.2638E − 5
3.1741
12.0453
1.7029E + 5
1.5373E + 5
5.6777E + 4
0.9751
112.2467
e
e
3.8143E + 6
3.7554E + 6
2.6576E + 6
2.0358E + 3
3.1395E + 3Implementation of DT in Matlab in case of Example2 produces numerical results which are listed in Table 3. The results of DT with order of the Taylor polynomial N=10 are compared to values obtained by DT combined with modified Adomian formula in [18] and to values produced by Matlab function DDENSD.Table 3
Example2, comparison of presented DT, DT in [18] and Matlab function DDENSD for u1 and u2 in the first step for t∈0,t1, and presented DT and DDENSD in the second step for t∈t1,t2.
Method
DT 10
DT 10
Matlab DDENSD
Matlab DDENSD
DT [18]
DT [18]
t
u
1
u
2
u
1
u
2
u
1
u
2
0.00
1.0000
0.0000
1.0000
0.0000
1.0000
0.0000
0.05
1.0513
0.0024
1.0513
0.0024
1.0513
0.0024
0.10
1.1054
0.0093
1.1050
0.0093
1.1051
0.0093
0.15
1.1625
0.0203
1.1614
0.0203
1.1618
0.0203
0.20
1.2230
0.0348
1.2204
0.0348
1.2209
0.0348
0.25
1.2871
0.0524
1.2822
0.0524
1.2832
0.0524
0.30
1.3552
0.0726
1.3469
0.0726
1.3481
0.0726
0.35
1.4277
0.0951
1.4146
0.0951
1.4160
0.0951
0.3520
1.4306
0.0960
1.4174
0.0960
—
—
0.3904
1.4896
0.1146
1.4717
0.1146
—
—
0.4289
1.5514
0.1340
1.5280
0.1340
—
—
0.4673
1.6161
0.1541
1.5864
0.1541
—
—
0.5057
1.6839
0.1747
1.6469
0.1748
—
—
0.5441
1.7548
0.1957
1.7098
0.1957
—
—
0.5826
1.8290
0.2169
1.7750
0.2169
—
—
0.6210
1.9065
0.2380
1.8428
0.2381
—
—First, we should say that the function DDENSD had difficulty at 0 where the value of the delayed argumentt/2 was equal to the argument itself. Hence, to make DDENSD work, we replaced t/2 by t/2−10−16 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 interval0,t2. When trying to follow the method of steps, i.e., using DDENSD on 0,t1 and then on t1,t2, the results on the second interval t1,t2 did not correspond to reality: there was a discontinuity in u2 at t1.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 componentu2, while we observe a growing distance between the values of the first component u1 computed by presented DT algorithm and values computed by the other two algorithms. As u1 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 0,t1 and t1,t2 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.


## 3.1. Algorithm Description

Recall system (1)(15)unt=ft,ut,u′t,…,un−1t,u1α1t,…,urαrt,with initial conditions(16)ut0=v0,u′t0=v1,…,un−1t0=vn−1,and initial vector function Φt on interval t∗,t0 satisfying(17)Φt0=ut0,…,Φn−1t0=un−1t0.Further recall that in Section2.2, we broke the interval t0,T into subintervals Il=tl−1,tl, l=1,…,K. Define I0=t∗,t0.Then, we are looking for a solutionut of the IVP (1)–(3) in the form(18)ut=uI1t,t∈I1,uI2t,t∈I2,⋮uIKt,t∈IK,where solution uIj in the jth interval Ij is obtained in the following way. We solve the following equation:(19)uIjnt=ft,uIjt,uIj′t,…,uIjn−1t,uj,1α1t,…,uj,rαrt,where(20)uj,iαit=uIlαit,uIl′αit,…,uIlmiαit,if
α
i
t
∈
I
l
,
for
t
∈
I
j
,
l
∈
1
,
…
,
j
,
j
∈
1
,
…
.
K
.In case thatαit∈I0=t∗,t0 for t∈Ij, then again(21)uj,iαit=ϕαit,ϕ′αit,…,ϕmiαit.Application of DT attj−1 to equation (19) yields a system of recurrence algebraic equations:(22)UIjk+ntj−1=Fk,UIjk,UIjk+1,…,UIjk+n−1,where the function F 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(23)UIjktj−1=1k!uIjktj−1,fork=0,1,…,n−1,j∈1,…,K.Using (22) with (23) and then inverse transformation rule, we obtain approximate solution to (19) in the form of Taylor series:(24)uIjt=∑k=0∞UIjktj−1t−tj−1k,t∈Ij,for all j∈1,…,K.To transform (20) correctly, we need the following theorem.Theorem 2.
Letαit∈Il for t∈Ij, where l∈1,…,j−1. Let p∈ℕ. Denote Dαittj−1=Aiktj−1k=0∞. Then,(25)DuIlpαittj−1=∑y=0ky+p!y!UIly+pαitj−1⋅B^k,yAi1tj−1,…,Aik−y+1tj−1k=0∞,where B^0,0=1, B^k,0=0 for k≥1, and(26)UIlyαitj−1=∑x=0∞x+yxαitj−1−tl−1xUIlx+ytl−1,for y≥0.Proof.
To prove (25) with p=0, we use Theorem 1 with ft=uIlt, gt=αit, and ht=f∘gt. We immediately get(27)Hktj−1=∑y=1kUIlyαitj−1⋅B^k,yAi1tj−1,…,Aik−y+1tj−1,for k≥1. For k=0, Theorem 1 yields H0tj−1=UIl0αitj−1=UIl0αitj−1⋅B^0,0Ai1tj−1. Now, (25) for p>0 is a consequence of Lemma 1 and it remains to prove (26). We recall that(28)uIlt=∑k=0∞UIlktl−1t−tl−1k,t∈Il.
As the assumption was thatαitj−1∈Il, we may apply Definition 2 to (28) and obtain(29)UIlyαitj−1=1y!dyuIltdtyt=t01y!∑z=y∞z!z−y!UIlztl−1t0−tl−1z−y.
Substitutingt0=αitj−1 and z=x+y gives (26).


## 3.2. New DT Formulas

In the applications, we also use the following DT formulas.Theorem 3.
Assume thatFkt0 is the differential transformation of the function ft and r∈R:(a)
Ifft=tr, then Fkt0=rkt0r−k for all t such that t−t0<t0,whererk=rr−1…r−k+1/k!=rk/k! and rk is the Pochhammer symbol.(b)
Ifft=lnt, then Fkt0=−1k−1/k⋅t0k for k≥1.Proof.
(a)
Recall the Newton’s generalisation of the binomial formula: ifx and y are real numbers with x>y, and r is any complex number, one has(30)x+yr=∑k=0∞rkxr−kyk,where rk=rr−1…r−k+1/k!. Let us rewrite tr as tr=t−t0+t0r=t0+t−t0r. Applying (30) yields(31)tr=∑k=0∞rkt0r−kt−t0k.(b)
We start by proving the formula(32)lntk=−1k−1k−1!tk,by induction. For k=1, we have lnt′=1/t; hence, (32) is valid. Suppose that (32) holds for k. Then,(33)lntk+1=lntk′=−1k−1k−1!tk′=−1k−1k−1!t−k′=−1k−1k−1!−kt−k−1=−1kk!tk+1.Thus, formula (32) is valid for all k∈ℕ. Now by Definition 2,(34)Fkt0=1k!dklntdtkt=t0=1k!−1k−1k−1!tkt=t0=−1k−1k⋅t0k.


## 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 Section3.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 functionet. 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:(35)u′t=ut−t+ulnt,with the initial condition(36)u1=e,and with the initial function(37)ϕt=et,t∈0,1.
First we find the differential transform of the initial condition (36) which is U01=e. Further denote Det0=Ek0k=0∞ as the transformation of the exponential function with the center at 0 and Dlnt1=Fk1k=0∞ as the transformation of the logarithmic function at 1, respectively. Then, Lemma 1 and Theorem 3 yield(38)Ek0=1k!,Fk1=−1k−1k,fork≥1.
Fort∈1,e, equation (35) is transformed into(39)k+1Uk+11=Uk1−δk−δk−1+Hk1,where(40)Hk1=∑l=1kEl0B^k,lF11,…,Fk−l+11,fork≥1,H01=E00=1.
We have(41)U11=e−1−0+1=e,H11=E10⋅F11=1⋅1=1,U21=12U11−0−1+H11=12e−1+1=e⋅12,H21=E10⋅F21+E20⋅F112=1⋅−12+12⋅1=0,U31=13U21+H21=13⋅e⋅12=e⋅13!,⋮
Using the inverse transformation, we see that fort∈1,e,(42)ut=e1+t−1+t−122+t−133!…=e⋅et−1=et,which corresponds to the exact solution to the IVPs (35)–(37).
In the second step of the method of steps, i.e., in the intervalt∈e,ee, we know that ut=et for t∈1,e and equation (35) is transformed into(43)k+1Uk+1e=Uke−eδk−δk−1+Hke,where(44)H0e=U01,Hke=∑l=1kUl1B^k,lF1e,…,Fk−l+1e,fork≥1.
Here,Fke, according to Theorem 3, are coefficients of Taylor series of logarithmic function with the center at e:(45)F0e=1,Fke=−1k−1k⋅ek,fork≥1.
Taking the values calculated in the first step and substituting them into the recurrence formulas (43) and (44), we obtain(46)U0e=ue=ee,U1e=U0e−eδ0−δ−1+H0e=ee−e−0+e=ee,H1e=U11⋅F1e=e⋅1e=1,U2e=12U1e−0−1+H1e=12ee−1+1=12ee,H2e=U11⋅B^2,1F1e,F2e+U21⋅B^2,2F1e=e⋅F2e+e2⋅F1e2=e⋅12−1e2+e2⋅1e2=0,U3e=13U2e+H2e=13!ee,H3e=U11⋅B^3,1F1e,F2e,F3e+U21⋅B^3,2F1e,F2e+U31⋅B^3,3F1e=e⋅F3e+e2⋅2⋅F1e⋅F2e+e3!F1e3=e⋅2e313!+e⋅1e⋅−1e2⋅12+e6⋅1e3=1e213−12+16=0,U4e=14U3e+H3e=14!ee.Hence, for t∈e,ee, we have(47)ut=ee+eet−e+12eet−e2+13!eet−e3+14!eet−e4+…=ee1+t−e+t−e22+t−e33!+t−e44!+…=ee⋅et−e=et,which 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:(48)u1‴=u1‴t−2u1t3+u1t23+u2′t−12e−t,u2‴=12u2‴t2+u2′t−1u1t3,with initial functions(49)ϕ1t=et,ϕ2t=t2,for t∈−2,0, and initial conditions(50)u10=1,u2′0=1,u1″0=1,u20=0,u2′0=0,u2″0=2.
Fort∈0,t1, where t1≈0,351734 is the minimal root of t−1/2e−t=0, using the method of steps, we obtain(51)u1‴=et−2u1t3+u1t23+2t−e−t,u2‴=12u2‴t2+2t−1u1t3.
We need to find the differential transform of the considered problem. We notice that system (2) contains nonlinear term ht=u1t23. To get DT of this term, Dht0=H1k0k=0∞, and we apply Theorem 1. First, applying DT to system (2) at t0=0, we get the recurrent system:(52)k+1k+2k+3U1k+30=e−2∑l=0k1l!13k−lU1k−l0+H1k0+2δk−−1kk!,(53)k+1k+2k+31−12k+1U2k+30=213k−1U1k−10−23kU1k0.
Denotegt=t2/3; then, ht=g∘u1t, and following Theorem 1, we obtain(54)H100=G101,H1k0=∑l=1kG1l1B^k,lU110,…,U1k−l+10,fork≥1, where Dgt1=G1k1k=0∞ and, Theorem 3 being applied, G1k1=2/3k for k≥0. Furthermore, the transformed initial conditions are(55)U100=1,U110=1,U120=12,U200=0,U210=0,U220=1.
Using them, we compute the first three coefficients of the nonlinear termht:(56)H100=G101=1,H110=G111⋅B1,1U110=23⋅1=23,H120=G111⋅B2,1U110,U120+G121⋅B2,2U110=23⋅12−19⋅1=29.
Solving recurrent systems (52) and (53), we get(57)k=0:U130=16e−2U100+H100+1=2+e−26,U230=13−2U100=−23.k=1:U140=e−22413U110+U100+124H110+1=4e−2+572,U240=1182U100−23U110=227.k=2:U150=e−26019U120+13U110+12U100+160H120−12=16e−2−51080,U250=210523U110−29U120=2189.
Using the inverse DT (Definition3), we get approximate solution for the IVPs (48)–(50) on the interval 0,t1:(58)u1,I1t=1+t+12t2+2+e−26t3+4e−2+572t4+16e−2−51080t5+…,u2,I1t=t2−23t3+227t4+2189t5+…,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 intervalt1,t2, where t2 is the minimal root of t−1/2e−t=t1, t2≈0,620556. Now taking into account that both proportional delays q1t=1/3t and q2t=1/2t and also the time-dependent delay t−τ1t=t−1/2e−t map the interval t1,t2 into the interval 0,t1, system (2) becomes(59)u1‴=et−2u1,I1t3+u1t23+u2,I1′t−12e−t,u2‴=12u2,I1‴t2+2t−1u1,I1t3.
DenoteDhtt1=H2kt1k=0∞ and Dftt1=F2kt1k=0∞, where ft=u2′t−1/2e−t. By application of Theorem 2 to corresponding terms, system (2) transformed at t0=t1 reads as(60)k+1k+2k+3U1k+3t1=e−2∑l=0ket1l!13k−l∑x=0∞x+k−lxU1x+k−l0t13x+H2kt1+F2kt1,(61)k+1k+2k+3U2k+3t1=12k+1k+3k+2k+1∑x=0∞x+k+3xU2x+k+30t12x+2∑l=0kt1−1δl+δl−113k−l∑x=0∞x+k−lxU1x+k−l0t13x.
Now denoteDgtu1t1=G2ku1t1k=0∞; then, according to Theorem 3, G2ku1t1=2/3ku1t12/3−k for k≥0 and Theorem 1 implies(62)H20t1=G20u1t1=u1t123,H2kt1=∑l=1kG2lu1t1⋅B^k,lU11t1,…,U1k−l+1t1,fork≥1.
Further denoteet=t−1/2e−t and Dett1=E2kt1k=0∞. Then, ft=u2′∘et and, since et1=0, Theorem 1 in combination with Lemma 1 yields(63)E2kt1=t1⋅δk+δk−1−12⋅e−t1−1kk!,k≥0,F20t1=U210=0,F2kt1=∑l=1kl+1U2l+10⋅B^k,lE21t1,…,E2k−l+1t1,fork≥1.
To get the initial dataU1kt1 and U2kt1 for k=0,1,2, we have to transform(64)uit=Ui00+Ui10t+Ui20t2+Ui30t3+…,at t1, i=1,2. For k=0,1,2, we have(65)Uikt1=Ui00+Ui101kt11−k+Ui202kt12−k+Ui303kt13−k+…,i.e.,(66)Ui0t1=∑k=0∞Uik0t1k,Ui1t1=∑k=0∞k+1Uik+10t1k,Ui2t1=∑k=0∞k+22Uik+20t1k.
The initial values att1 will be approximated by taking finite sums in computer evaluations of the infinite sums above. Observe that u1t1=U10t1.
Now let us compute the first few values ofH2kt1. Denote B^k,l=B^k,lU11t1,…,U1k−l+1t1. Then, the first values of B^k,l are(67)B^0,0=1,B^1,1=U11t1⋅B^0,0=U11t1,B^2,1=12U11t1⋅B^1,0+U12t1⋅B^0,0=U12t1,B^2,2=U11t1⋅B^1,1=U11t12,and the coefficients H2 for k=0,1,2 are(68)H20t1=u1t1,H21t1=23u1t1−1/3⋅B^1,1=23u1t1−1/3⋅U11t1,H22t1=23u1t1−1/3⋅B^2,1−19u1t1−4/3⋅B^2,2=23u1t1−1/3⋅U12t1−19u1t1−4/3⋅U11t12.
Let us turn our attention to the first few values ofF2kt1. Starting with E2kt1,(69)E20t1=t1−12e−t1=0,E21t1=1+12e−t1=1+t1,E22t1=−12⋅e−t12=−t12.
Now,B^k,l are B^k,lE21t1,…,E2k−l+1t1:(70)B^0,0=1,B^1,1=E21⋅B^0,0=1+t1,B^2,1=12E21⋅B^1,0+E22⋅B^0,0=−t12,B^2,2=E21⋅B^1,1=1+t12.
Finally, coefficientsF2 for k=0,1,2 are(71)F20t1=U210=0,F21t1=2U220⋅B^1,1=21+t1,F22t1=2U220⋅B^2,1+3U230⋅B^2,2=−t1−21+t12.
At this moment, we substituteH2 and F2 into systems (60) and (61). The next three coefficients at t1 for U1 are(72)k=0:U13t1=16e−2et1⋅u1t13+u1t1,k=1:U14t1=124e−2et1⋅13⋅∑x=0∞x+1U1x+10t13x+et1⋅u1t13+23⋅u1t1−1/3⋅U11t1+21+t1,k=2:U15t1=160e−2et1132⋅∑x=0∞x+2xU1x+20t13x+et1⋅13⋅∑x=0∞x+1U1x+10t13x+et12⋅u1t13+16023u1t1−1/3⋅U12t1−19u1t1−4/3⋅U11t12−t1−21+t12,and for U2, we obtain(73)k=0:U23t1=12∑x=0∞x+3xU2x+30t12x+13t1−1⋅u1t13,k=1:U24t1=14∑x=0∞x+4xU2x+40t12x+112t1−1⋅13⋅∑x=0∞x+1U1x+10t13x+u1t13,k=2:U15t1=18∑x=0∞x+5xU2x+50t12x+130t1−1⋅19⋅∑x=0∞x+2xU1x+20t13x+13∑x=0∞x+1U1x+10t13x⋅
Using the inverse DT, again we get approximate solution for the IVPs (48), (49), and (50) on the intervalt1,t2:(74)u1,I2t=U10t1+U11t1t−t1+U12t1t−t12+U13t1t−t13+…,u2,I2t=U20t1+U21t1t−t1+U22t1t−t12+U23t1t−t13+….
As the calculations are getting more complicated, all the calculations have been done numerically only.


## 3.4. Numerical Results and Discussion

Table1 shows comparison of results for Example 1 obtained by DT algorithm with the orders of Taylor polynomials of the approximate solution N=5,10,25 to results of Matlab function DDENSD in the interval 1,e. 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 N=10,25 does even better, whereas DT for N=5 does not show satisfactory precision.Table 1
Example1, error analysis of u in 1,e.
Exact solution
DT 5
DT 10
DT 25
Matlab DDENSD
t
e
t
u
−
e
t
u
−
e
t
u
−
e
t
u
−
e
t
1
2.7183
0
0
0
0
1.4296
4.1769
0.0000
6.4890E − 12
8.8818E − 16
1.6976E − 5
1.8591
6.4182
0.0017
1.0381E − 8
8.8818E − 16
3.6095E − 5
2.2887
9.8622
0.0211
1.2409E − 6
3.5527E − 15
9.3438E − 5
e
15.1543
0.1273
3.0578E − 5
1.2861E − 12
1.4012E − 4Table2 brings the same comparison in the second interval e,ee. 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 e,ee, the rate of precision would be better seen using relative errors.Table 2
Example1, error analysis of u in e,ee.
Exact solution
DT 5
DT 10
DT 25
Matlab DDENSD
t
e
t
u
−
e
t
u
−
e
t
u
−
e
t
u
−
e
t
e
15.1543
3.5527E − 15
3.5527E − 15
3.5527E − 15
0
5.8273
339.4331
32.3491
1.333
2.2737E − 13
0.0743
8.9363
7.6028E + 3
4.4755E + 3
397.5509
2.2638E − 5
3.1741
12.0453
1.7029E + 5
1.5373E + 5
5.6777E + 4
0.9751
112.2467
e
e
3.8143E + 6
3.7554E + 6
2.6576E + 6
2.0358E + 3
3.1395E + 3Implementation of DT in Matlab in case of Example2 produces numerical results which are listed in Table 3. The results of DT with order of the Taylor polynomial N=10 are compared to values obtained by DT combined with modified Adomian formula in [18] and to values produced by Matlab function DDENSD.Table 3
Example2, comparison of presented DT, DT in [18] and Matlab function DDENSD for u1 and u2 in the first step for t∈0,t1, and presented DT and DDENSD in the second step for t∈t1,t2.
Method
DT 10
DT 10
Matlab DDENSD
Matlab DDENSD
DT [18]
DT [18]
t
u
1
u
2
u
1
u
2
u
1
u
2
0.00
1.0000
0.0000
1.0000
0.0000
1.0000
0.0000
0.05
1.0513
0.0024
1.0513
0.0024
1.0513
0.0024
0.10
1.1054
0.0093
1.1050
0.0093
1.1051
0.0093
0.15
1.1625
0.0203
1.1614
0.0203
1.1618
0.0203
0.20
1.2230
0.0348
1.2204
0.0348
1.2209
0.0348
0.25
1.2871
0.0524
1.2822
0.0524
1.2832
0.0524
0.30
1.3552
0.0726
1.3469
0.0726
1.3481
0.0726
0.35
1.4277
0.0951
1.4146
0.0951
1.4160
0.0951
0.3520
1.4306
0.0960
1.4174
0.0960
—
—
0.3904
1.4896
0.1146
1.4717
0.1146
—
—
0.4289
1.5514
0.1340
1.5280
0.1340
—
—
0.4673
1.6161
0.1541
1.5864
0.1541
—
—
0.5057
1.6839
0.1747
1.6469
0.1748
—
—
0.5441
1.7548
0.1957
1.7098
0.1957
—
—
0.5826
1.8290
0.2169
1.7750
0.2169
—
—
0.6210
1.9065
0.2380
1.8428
0.2381
—
—First, we should say that the function DDENSD had difficulty at 0 where the value of the delayed argumentt/2 was equal to the argument itself. Hence, to make DDENSD work, we replaced t/2 by t/2−10−16 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 interval0,t2. When trying to follow the method of steps, i.e., using DDENSD on 0,t1 and then on t1,t2, the results on the second interval t1,t2 did not correspond to reality: there was a discontinuity in u2 at t1.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 componentu2, while we observe a growing distance between the values of the first component u1 computed by presented DT algorithm and values computed by the other two algorithms. As u1 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 0,t1 and t1,t2 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.Example1 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 Example2, we could confirm a very good concordance of both differential transform and DDENSD in values of the component u2 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 u1 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.



---

*Source: 2854574-2020-02-28.xml*