A Numerical Method for Delayed Fractional-order Differential Equations

A numerical method for nonlinear fractional-order differential equations with constant or time-varying delay is devised. The order here is an arbitrary positive real number, and the differential operator is with the Caputo definition. The general Adams-Bashforth-Moulton method combined with the linear interpolation method is employed to approximate the delayed fractional-order differential equations. Meanwhile, the detailed error analysis for this algorithm is given. In order to compare with the exact analytical solution, a numerical example is provided to illustrate the effectiveness of the proposed method.


Introduction
Fractional calculus is an old mathematical problem and has been always thought of as a pure mathematical problem for nearly three centuries [1][2][3]. Though having a long history, it was not applied to physics and engineering for a long period of time. However, in the last few decades, fractional calculus began to attract increasing attention of scientists from an application point of view [4][5][6]. In the fields of continuous-time modeling, many scholars have pointed out that fractional derivatives are very helpful in describing linear viscoelasticity, acoustics, rheology, polymeric chemistry, and so forth. Moreover, fractional derivatives have proved to be a very suitable tool for the description of memory and hereditary properties of various materials and processes. Nowadays the mathematical theories and practical applications of these operators are well established, and their applicabilities to science and engineering are being considered as attractive topics. In modeling dynamical processes, equations described with fractional derivatives have appeared in various fields such as physics, chemistry, and engineering [7,8], and therefore, fractional-order differential equations (FDEs) are becoming the focus of many studies.
The development of effective methods for numerically solving FDEs has received increasing attention over the last few years. Several methods based on Caputo or Riemann-Liouville definition have been proposed and analyzed [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. For instance, based on the predictor-corrector scheme, Diethelm et al. introduced Adams-Bashforth-Moulton algorithm [11][12][13], and meanwhile some error analysis and an extension of Richardson extrapolation were presented to improve the numerical accuracy. These techniques elaborated in [11][12][13] take advantage of the fact that the FDEs can be reduced to Volterra type integral equations. And therefore, one can apply the numerical schemes for Volterra type integral equations to find the solution of FDEs; one can refer to [21,[25][26][27][28] and the literature cited therein for more details. The technique presented in [11][12][13] has been further analyzed and extended for multiterm [15] and multiorder systems [17]. Li and Tao [20] studied the error analyses of the fractional Adams method for the fractional-order ordinary differential equations for more general cases. Deng [24] obtained a good numerical approximation by combining the short memory principle with the predictor-corrector approach.
However, in real world systems, delay is very often encountered in many practical systems, such as automatic control, biology and hydraulic networks, economics, and long transmission lines. Delayed differential equations are correspondingly used to describe such dynamical systems. In recent years, delayed FDEs begin to arouse the attention 2 Journal of Applied Mathematics of many researchers [29][30][31]. To simulate these equations is an important technique in the research, accordingly, finding effective numerical methods for the delayed FDEs is a necessary process. However, to the best of our knowledge, there are very few works devoted in the literature to such work so far. The aim of this paper is to generalize the Adams-Bashforth-Moulton algorithm introduced in [11][12][13] to the delayed FDEs. The rest of this paper is organized as follows. In Section 2, the numerical scheme is devised based on the Adams-Bashforth-Moulton method, and meanwhile the predictorcorrector scheme is developed in the case of constant and time-varying delay, respectively. The detailed error analysis of the newly proposed algorithm is given in Section 3. Numerical simulations are performed in Section 4 in order to illustrate the effectiveness of the proposed scheme. Finally, some concluding remarks are reported in Section 5.

Formulation of the Numerical Method for Delayed FDEs
There are three commonly used definitions for fractional derivatives [2], that is, Riemann-Liouville, Grünwald-Letnicov, and Caputo definitions. In pure mathematics, Riemann-Liouville derivative is more commonly used than Caputo derivative. However, as we all know, only the Caputo definition has the same form as integer-order differential equations in the initial conditions. And therefore this paper is based on Caputo definition.

Caputo Fractional Derivatives
Definition 1 (see [2]). The Caputo fractional derivative of order of a continuous function : + → is defined as follows: where Γ is Γ-function, and

The Procedure of Adams-Bashforth-Moulton Algorithm.
The main aim of this paper is to study a numerical scheme for the approximate solution of delayed FDEs. For this purpose, we consider delayed FDEs described by where is the order of the differential equation, ( ) is the initial value, and is an integer. Before introducing the new algorithm for the delayed FDEs, we first briefly recall the idea of the classical one-step Adams-Bashforth-Moulton algorithm for the following ODE: Let ∈ [0, ], = ℎ, = 0, 1, . . . , , and ℎ = / be the time step. Assume that = ( ), = 1, . . . , , then the previous equation is equivalent to Since the integral in the right-hand side of (6) can be approximated by trapezoidal quadrature formula, ∫ +1 ( , ( )) ≈ ℎ 2 ( ( , ( )) + ( +1 , ( +1 ))) .
Thus, we can get the approximation to +1 as follows: which is an implicit equation and cannot be solved for +1 directly; however, one can adopt another numerical method to approximate +1 in the right-hand side preliminarily. The preliminary approximation +1 is called predictor and can be obtained by forward Euler method or any other explicit method; here we take forward Euler method as an example; that is, Then, the so-called one-step Adams-Bashforth-Moulton method is  in which (9) is called the predictor, and (10) is called the corrector. It is well known that this method is convergent; that is, max =0,1,...,

The Representation of Algorithm for Delayed FDEs.
Having introduced the Adams-Bashforth-Moulton algorithm, now the key problem is to establish approximation to the delayed term ( − ), which consists of the following two cases.
Case I (when is constant). It can be seen that, for any positive constant , − may not be a grid point for any . Suppose that ( − )ℎ = and 0 ≤ < 1. When = 0, ( − ) can be approximated by and when 0 < < 1, ( − ) cannot be calculated directly. For ( − 1)ℎ < < ℎ, let V +1 be the approximation to ( +1 − ); we interpolate it by the two nearest points; that is, It can be seen from (13) that if > 1, the numerical equation is explicit and thus can be computed directly. However, when = 1 and ̸ = 0, that is, < ℎ, then the first term in the righthand side of the previous equation is +1 . It still needs to predict, in this case, that is Case II (when is time varying). If is time varying, that is, = ( ), then the approximation in this case is more complicated. Let V +1 be the approximation to ( +1 − ), and the linear interpolation of { } at point = +1 − ( +1 ) is employed to approximate the delay term. Let ( +1 ) = ( +1 − +1 )ℎ, where +1 is a positive integer and +1 ∈ [0, 1), then It should be noted that when is constant, for given ℎ and , we can judge whether = 1 or > 1 holds at the beginning of the programme. However, when is time varying, is also time varying; that is, at one moment it is equal to 1, and at another moment it may greater than 1. If +1 = 1, the first term in the right-hand side of (15) still needs to predict; otherwise, when +1 > 1, it needs not. So in each step of the computing, it needs to test whether = 1 or not firstly, then computation of the first term in the right-hand side of (15) depends on whether it requires predicting or not in different steps.

Error Analysis of the Newly-Proposed Algorithm
Before giving the main results, we firstly introduce the following two lemmas, which will be used in the proof of the theorem.
Proof. We will prove the result by virtue of mathematical induction. Suppose that the conclusion is true for = 0, 1, . . . , , then we have to prove that the inequality also holds for = + 1. To do this, we first consider the error of the predictor +1 .
From the assumptions (24) So and from we have Since we have which completes the proof.
When 1 = 1, 2 = 2, applying Theorem 3, it is not difficult to derive the following result on the error bound. Remark 5. We can see from Cases I and II that if the solution ( ) is not smooth, then the truncation error of V +1 , which approximates to ( − ), will not be (ℎ 2 ). Therefore, by the error analysis elaborated in Section 3, we cannot yield the same conclusion as in Theorem 4. We will consider such a problem in our future work.

A Numerical Example
In this section, the following delayed FDE is considered: Notice that the exact solution to this equation is In accordance with delay being constant or time varying, the numerical results are displayed in Tables 1, 2, 3, and  4. respectively, where denotes the absolute numerical error and denotes the relative numerical error. From the numerical results we can see that the computing errors are in general acceptable for engineering.

Conclusions
In this paper, a numerical algorithm is formulated based on Adams-Bashforth-Moulton method for fractional-order differential equations with time delay. The error analysis of the numerical scheme is carried out, meanwhile, a numerical example with constant delay and time-varying delay is proposed to testify the effectiveness of the scheme. This algorithm can be used not only in the simulation of delayed fractional-order differential equations but also in the simulation of delayed fractional-order control systems.