Abstract

Fractional functional differential equations with delay (FDDEs) have recently played a significant role in modeling of many real areas of sciences such as physics, engineering, biology, medicine, and economics. FDDEs often cannot be solved analytically so the approximate and numerical methods should be adapted to solve these types of equations. In this paper we consider a new method of backward differentiation formula- (BDF-) type for solving FDDEs. This approach is based on the interval approximation of the true solution using the Clenshaw and Curtis formula that is based on the truncated shifted Chebyshev polynomials. It is shown that the new approach can be reformulated in an equivalent way as a Runge-Kutta method and the Butcher tableau of this method is given. Estimation of local and global truncating errors is deduced and this leads to the proof of the convergence for the proposed method. Illustrative examples of FDDEs are included to demonstrate the validity and applicability of the proposed approach.

1. Introduction

Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences [14]. Some approximate and numerical methods are used to obtain approximate solutions [511]. The general theory of differential equations with delays (DDEs) is widely developed and discussed in the literature [1216]. DDEs are differential equations in which the rate of change of does not depend only on the values of for the same time value but also on time values less than . In the simplest case, DDEs have the formunder the following initial conditions: where is the delay term. DDEs arise in scientific models. These models of DDEs appeared in [13, 14, 17, 18]. Recently, fractional differential equations with delay (FDDEs) gain the attention of many researchers. FDDEs have appeared also in some scientific areas [1923]. Some approaches are built to obtain numerical and approximate solutions for FDDEs:under the following initial conditions: Yang and Cao [24] studied the existence and uniqueness of initial value problems for nonlinear higher fractional equations with delay by fixed point theory. The theory of fractional functional differential equation and the asymptotic properties of fractional delay differential equations are discussed, respectively, in [2527]. Bhalekar and Daftardar-Gejji adapted the Adams Bashforth Moulton predictor corrector method (EABMPC method) [28] to solve this form of FDDEs (3)-(4), where . Besides, Wang [29] combined Adams Bashforth Moulton method with the linear interpolation method to approximate FDDEs. In addition, Wang et al. [30] introduced a numerical method based on Grunwald Letnikov definition to solve nonlinear FDDEs with constant time varying delay.

In addition to that, Morgado et al. [31] made an analysis and numerical methods for linear fractional differential equations with positive finite delay: under the following initial conditions: Also, Moghaddam and Mostaghim [32] discussed and introduced a novel matrix approach to fractional finite difference for solving models based on nonlinear fractional delay differential equations of the form (5)-(6). Also, Moghaddam and Mostaghim [33] developed a numerical method based on finite difference for solving fractional delay differential equations of the form (5)-(6). In [34] the authors consider a fourth-order method of BDF-type for solving stiff initial-value problems, based on the interval approximation of the true solution by truncated Chebyshev series. It is shown that the method may be formulated as a Runge-Kutta method having stage order four.

In this paper, we focus on FDDEs with finite delay which has the following form:We present fractional order Runge-Kutta method of backward differentiation formula- (BDF-) type based on approximations by truncated shifted Chebyshev series of Clenshaw and Curtis formula in the sense that we do not take an approximation of the function in the right-hand side of the differential equation, which will be integrated later, but we take a truncated Chebyshev series instead.

The structure of this paper is arranged in the following way: in Section 2, we present some necessary definitions and mathematical preliminaries of the fractional calculus theory. In Section 3, we introduce the derivation of the method and the fundamental related relations. In Section 4, we present the proposed method as a one-step recurrence formula, obtain the local truncating error, and prove that the order of the proposed method is such that is the approximation order of Clenshaw and Curtis formula. In Section 5, the global truncating error of our proposed method is deduced. In Section 6, numerical examples are given to solve FDDEs and show the accuracy of the presented method. Finally, in Section 7, the paper ends with a brief conclusion and some remarks.

2. Preliminaries and Notations

In this section, we present some necessary definitions and mathematical preliminaries of the fractional calculus theory required for our subsequent development.

2.1. The Caputo Fractional Derivatives

Definition 1. The Caputo fractional derivative operator of order is defined in the following form: where , , .

Similar to integer-order differentiation, the Caputo fractional derivative operator is a linear operation: where and are constants.

For Caputo’s derivative we have We use the ceiling function to denote the smallest integer greater than or equal to and . Recall that, for , the Caputo differential operator coincides with the usual differential operator of integer order.

For more details on fractional derivatives definitions and their properties, see [3538].

2.2. The Definition and Properties of the Shifted Chebyshev Polynomials

The well-known Chebyshev polynomials [39] are defined on the interval and can be determined with the aid of the following recurrence formula:It is well known that , . The analytic form of the Chebyshev polynomials of degree is given by where denotes the integer part of . The orthogonality condition isIn order to use these polynomials on the interval we define the so-called shifted Chebyshev polynomials by introducing the change of variable . The shifted Chebyshev polynomials are defined asThe analytic form of the shifted Chebyshev polynomial of degree is given by where and . The orthogonality condition of these polynomials is where the weight function , , with , , .

The function , which belongs to the space of square integrable in , may be expressed in terms of shifted Chebyshev polynomials as where the coefficients are given by The well-known shifted Chebyshev polynomials of the first kind of degree are defined on the interval as in (15). We choose the grid (interpolation) points to be the Chebyshev-Gauss Lobatto points associated with the interval , , . These grids can be written as .

Clenshaw and Curtis [40] introduced an approximation of the function , and we reformulate it to be used on the shifted Chebyshev polynomials as follows: The summation symbol with double primes denotes a sum with both first and last terms halved.

3. Derivation of the Numerical Approximation Scheme

Let us consider the initial value problem for arbitrary order fractional differential equation with delay of the formsuch that and is an integer.

Yang and Cao studied the existence and uniqueness of these types of (20) using the fixed point theory [24]. Now, let be an approximation of the theoretical solution at . We are interested in obtaining a numerical approximation at the point such that is the step size.

If we rewrite the solution in the interval in terms of a new variable defined bywe can approximate the solution of FDDEs expressed in the form Approximate (22) in terms of Clenshaw and Curtis finite sum approximations as follows: such that Using the analytic form of shifted Chebyshev polynomials (15) and the properties of Caputo fractional derivatives (10) givesNow, can be expressed approximately in terms of shifted Chebyshev series as follows: where is obtained from (18) as clarified by Doha et al. [5],From (27)–(29), we getwhere ,Then, can be put in the following form:Then,From (23) and (33), we obtain the following approximate equation: Evaluating the formula (34) at , gives an implicit system of algebraic equations

It can be observed that, for any delayed term , may not be a grid point for any . So, we establish the approximation for the delayed function as follows:So, we define the delayed term as such that and is a positive integer. As , then Then, we have two different cases.

The First Case (). Considersuch that

The Second Case (). In this case, is not a grid point. To overcome this, we introduce the following approximation:Then, we have such that From (35)–(42), we have the following numerical approximation scheme:where the unknowns are the values of the solution at the intermediate points .

For DDEs , explicitly, if we put in (43), the algebraic system will have the following form:For FDDEs , also, if we put in (43), the algebraic system will have the following form:such thatFor the sake of simplicity, we take this abbreviation Solving this system, we obtain in particular the required value at the final point on the interval . If we repeat this procedure along the integration interval , a discrete solution for the problem in (20) will be deduced.

4. One-Step Recurrence Formula Expression and the Local Truncating Error

The method in (43) at can be expressed by one step recurrence formula: where

Remark 2. ConsiderFor DDEs , the constant matrices and are given by which have a completely agreement with the work done by Ramos and Vigo-Aguiar [34].

For FDDEs , the constant matrices and are also given byAs the matrix is nonsingular, multiplying the formula (48) by the inverse matrix yields

Remark 3. The product is a square matrix of dimension on the following form:

Remark 4. The method expressed finally by (53) can be considered as a Runge-Kutta method where the Butcher tableau [41] is given by the following.(i)For DDEs , see Table 1.When we compare this tableau with the tableau obtained by Ramos and Vigo-Aguiar [34], we find a complete agreement.(ii)For FDDEs , see Table 2.The contents of Butcher tableau at can be written also in terms of fractions but it will be complicated. So, we write it in decimal form.
In the names of Remark 4, we can say that the theory of Runge-Kutta methods may be applied here. For our approach, when we deal with DDEs , the Butcher conditions for order four are verified. So, this method has the fourth order for DDEs.

Consider the difference operators associated with each stage in our method as follows: and using the fractional Taylor series about , we obtain

Theorem 5. If the proposed Runge-Kutta method with shifted Chebyshev polynomials (48) is of order and if all partial fractional derivatives of up to order exist and are continuous, then the local error of (56) admits the rigorous bound:and so

Proof. Because of the order conditions and using (57), the first terms in the expansion of (56) vanish. Then (58) is obtained and the proof is complete.

Under the assumption of Theorem 5, we see that the expressions (58) and (59) are bounded by a constant independent of .

Remark 6. The proposed Runge-Kutta method with shifted Chebyshev polynomials (48) has order for sufficiently smooth problem FDDEs (20) if (59) is achieved.

5. Estimation of the Global Truncating Error

The global error is the error of the approximated solution after many steps.

Theorem 7. Suppose that is the exact solution of the system of FDDEs (20) and that is the approximate solution. If the following inequalities are held:then for , one has the following error estimate:

Proof. For any chosen norm, we investigate the error estimated by and now, we try to estimate its natural growth.
Operate on both sides of (64) by the Caputo fractional derivative operator and using the triangle inequality yieldsfrom Lipchitz condition (62), and we introduce which is called the defect of the approximate solution such that ; we have Now, we are about to solve instead of (67) the fractional differential equation and conclude that To solve (69)-(70), we use the fractional complex transform and the modified Riemann-Liouville fractional derivative defined bywhere and .
Modified Riemann-Liouville fractional derivative has the following property: Then (69)-(70) are transformed toThe exact solution for (73) can be written as Then the solution can be written using the original variables as follows:A combination of (68) and (75) yields the desired result (63).

Theorem 8. Let be a neighborhood of where is the exact solution of (20). Suppose that in and that the local error estimates (59) are valid in . Then the global error can be estimated by where ,And is small enough for the numerical solution to remain in .

Proof. Our task is to estimate the global error and this estimate is found along steps of the numerical method (see Figure 1).
Estimate the local error using Theorem 5 and Figure 1 to obtain From Theorem 7 with , we obtain from Figure 1, ,as , then we can deduce that and from (81)–(83), we deduce thatThenA combination of (84)-(85) yields the desired result (77).

Remark 9 (see [43]). Since the global truncating error (77) tends to zero as tends to zero, then the proposed method (43) is convergent; that is, as .

6. Applications and Numerical Simulation

In this section, we illustrate the theoretical results obtained in the previous sections. In order to do this we will consider the following examples whose exact solutions are known.

Example 1. Consider The exact solution for this problem is .

For FDDEs (86) at (, , ) and (, , ) the absolute errors between the computed points and exact ones are presented in Figures 2 and 3 respectively.

Example 2. Consider where is defined by such that is a special case of hypergeometric function.
The exact solution for this problem is .
In Table 3, , we present the results obtained with the proposed method at for different step sizes such that the maximum of the absolute error is given by In the third column, the numerical convergence order which is defined by is introduced.

Example 3. Consider where is defined byThe exact solution for this problem is .
In Table 4, the computational results of the maximum of absolute errors are displayed at for different time steps with . In the third column, the numerical convergence order is introduced also. These results are in concordance with our theoretical results.

7. Conclusion and Remarks

In this paper, a new fractional order Runge-Kutta method based on BDF-type Chebyshev approximations is introduced. This approach is applicable to initial value problems for arbitrary order fractional differential equations with delay. This new method can be expressed as one-step recurrence formula. It is shown that the method may be formulated in an equivalent way as a Runge-Kutta method of order . The local and global truncating errors for this new scheme are obtained and proved. Numerical examples with constant delay or time varying delay are proposed to justify the effectiveness of the proposed scheme and find a good agreement with the theoretical results. For the future work, the method is also applicable to systems of fractional functional differential equations with delay.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by Act 211 Government of the Russian Federation Program 02.A03.21.0006 on 27.08.2013 and by RFBR Grant 13-01-00089.