An Efficient Numerical Scheme for Solving a Fractional-Order System of Delay Differential Equations.

Fractional order systems of delay differential equations are very advantageous in analyzing the dynamics of various fields such as population dynamics, neural networking, ecology, and physiology. The aim of this paper is to present an implicit numerical scheme along with its error analysis to solve a fractional-order system of delay differential equations. The proposed method is an extension of the L1 numerical scheme and has the error estimate of O ( h 2 ) , where h denotes the step size. Further, we solve various non-trivial examples using the proposed method and compare the results with those obtained by some other established methods such as the fractional Adams method and the three-term new predictor-corrector method. We observe that the proposed method is more accurate as compared to the fractional Adams method and the new predictor-corrector method. Moreover, it converges for very small values of the order of fractional derivative.

Introduction

Fractional Calculus (FC), where the derivatives and integrals are considered in arbitrary form is one of the most currently active areas of research. The genesis of FC began with a question raised by L’ Hôpital to Leibniz towards the end of the 17th century. However, in the midway of the 19th century, the pioneering works of various mathematicians such as Liouville, Riemann, Grunwald, and Letnikov led to the formulation of fractional integrals and derivatives with the subsequent development of FC [1, 2]. Fractional derivatives are non-local in nature and preserved demonstrative physical characteristics. Hence, it is easy to anticipate and evaluate the dynamics of the various natural phenomena. However, FC has a long history but, its utility and capabilities have been realized in the past two-three decades. The applications of FC has been noticed in diverse fields such as earth system dynamics, epidemiology, computer vision, biot-theory, robotics, soil hydrology and mechanics, criminology, artificial neural networks, thermodynamics [3-12], and many other branches of science and engineering [13, 14]. Fractional delay differential equations (FDDEs) are the most suitable tools to model various real-world problems that incorporate history. Many authors have discussed the existence and uniqueness of FDDEs in the literature [15-17]. Moreover, widespread applications of FDDEs have been realized in various fields such as infectious diseases, immune systems, epidemics, tumor growth, population dynamics, circulating blood the body’s reaction to carbon dioxide, ecology, physiology [18]. Recently, various models based on FDDEs such as corona-virus disease model [19], hand-foot-mouth disease model [20], glucose-insulin interaction model [21] and so on have been discussed in the literature. FDDEs are more complex due to the involvement of fractional derivatives and delay terms. Thus, it is a challenging task to solve these equations analytically. Moreover, FDDEs do not possess exact solutions in most cases. Hence, one must depend on numerical methods. Therefore, the development of new accurate numerical methods for solving FDDEs is highly necessitated. In this pursuit, various classical methods have been modified and extended to solve FDDEs such as the fractional Adams method (FAM) [22], decomposition, and iterative methods [23, 24], operational matrix-based method [25], Runge–Kutta methods [26], wavelet methods [27, 28], Adams–Bashforth–Moulton method [29], finite difference methods [30], new predictor–corrector method (NPCM) [31], shifted Jacobi polynomial method [32], and so on. One of the finite difference methods is L1, where the fractional derivatives are discretized. The method L1 has been used by many researchers directly or indirectly to solve differential equations of fractional orders [33, 34]. In this paper, we extend the L1 numerical method for solving FDDEs of the following form:where is a constant delay and the operator denotes Caputo derivative defined as follows:Further, we present its error analysis. Furthermore, we exhibit the utility and applicability of the proposed method by performing some of the numerical simulations corresponding to chaotic and non-chaotic systems of FDDEs. We compare the solutions with exact, FAM, and NPCM. We observe that the proposed method has higher accuracy than FAM and NPCM. Moreover, the proposed method converges for very small values of , when FAM and NPCM both fail to converge.

This paper is organized as follows: In Sect. 2, we extend the L1 numerical method for solving FDDEs and present its error analysis in Sect. 3. In Sect. 4, we use the proposed method to solve various non-linear systems of FDDEs and compare the results with FAM and NPCM. Finally, in Sect. 5, we draw the conclusions.

Formulation of Numerical Method

In this section, we develop a numerical method for solving fractional delay differential equations (FDDEs). Consider the following general form of a non-linear system of FDDEs:where and are known functions and a delay term.

Consider a partition of the interval [0, T] with uniform grid points where M and K are positive integers such that and . Further, let and .

A numerical scheme to solve the system (1) is devised as follows: The Caputo derivative at is discretized as follows (cf. L1 algorithm [35]):whereSuppose represents the approximate value of at . The following process is used to compute the n-th approximation , while for are already computed. Approximate the fractional Caputo-derivative term that appear in Eq. (1) by the expression (2), we getwhere denotes the approximate value of at . After simplifying Eq. (4), we obtainFurther simplifying Eq. (5), we getOrAfter simplifying Eq. (7), we obtainSet:Thus, in view of the Eqs. (8) and (9), we obtain the following numerical scheme to solve the system of FDDEs (1):where are calculated by the formulas given in Eq.’s (3) and (9).

Error Analysis

In this section, we establish the error analysis of the proposed numerical method.

Lemma 1

[33, 36, 37] For , the following inequality holdswhere and is a constant.

We set:In view of the Eq.’s (11) and (12), we get

Theorem 1

Suppose satisfy the Lipschitz condition such that , where and are the Lipschitz constants. Let be the exact solution of the system (1) and the approximate solution at obtained by the proposed numerical method (10). Subsequently, for a sufficiently small value of h, we have

Proof

The numerical scheme given in Eq. (10) can be written aswhich is equivalent toSuppose that for Therefore, Eq. (15) can be written as:On account of Eq.’s (10) and (12), we haveIn view of Eq. (17), Eq. (16) turns out to beSince , therefore Eq. (18) implieswhere C is an arbitrary constant and does not depend on h. After simplifying Eq. (19), we obtainAs h is sufficiently small, we havewhich is the desired result.

Comment: We proved that the error estimate of the proposed method is . Whereas, the error-estimate of FAM is and for NPCM it lies between and . Hence, the proposed method has a better error estimate as compared to FAM and NPCM.

Illustrative Examples

In this section, we demonstrate the applicability of the proposed method by solving some of the non-trivial systems of FDDEs.

Example 1

Consider the following fractional delay differential equation:

where . Exact solution of the system (20) is . For and , we solve the system (20) using FAM, NPCM and the proposed method (10). Absolute and relative errors are calculated at and compared with FAM and NPCM in Tables 1, 2 and Figs. 1, 2. Obtained results reveal that the proposed method is more accurate than FAM and NPCM. Further, we compute the execution time taken by FAM, NPCM, and the present method. We found that FAM, NPCM, and the present method take 466.51, 234.57, and 425.93 s respectively. Hence the present method is more time-efficient as compared to FAM.

Table 1

Comparison of absolute errors for the system (20)

t	FAM	NPCM	Present method (10)
1	0.004391972	0.004392520	0.004117906
2	0.004115289	0.004116835	0.003470192
3	0.003961597	0.003964149	0.002966255
4	0.003855740	0.003859303	0.002522327
5	0.003775297	0.003779873	0.002112362
6	0.003710572	0.003716162	0.001724671
7	0.003656503	0.003663109	0.001352938
8	0.003610120	0.003617743	0.000993328
9	0.003569535	0.003578177	0.000643324
10	0.003533484	0.003543134	0.000301175

Table 2

Comparison of relative errors for the system (20)

t	FAM	NPCM	Present method (10)
2	0.002057644	0.002058418	0.001735096
3	0.000660266	0.000660691	0.000494376
4	0.000321312	0.000321608	0.000210194
5	0.000188765	0.000188994	0.000105618
6	0.000123686	0.0001238720	5.7489×10^-5
7	8.7059×10^-5	8.7217×10^-5	3.2213×10^-5
8	6.4466×10^-5	6.4603×10^-5	1.7748×10^-5
9	4.9577×10^-5	4.9697×10^-5	8.9351×10^-6
10	3.9260×10^-5	3.9368×10^-5	3.3464×10^-6

Fig. 1

Absolute errors (20)

Fig. 2

Relative errors (20)

Example 2

Consider the following fractional delay differential equation:

Exact solution of the system (21) is . We solve the system (21) for very small values of i.e. for and by using the proposed method, FAM and NPCM. At the numerical solutions obtained by these methods are compared in Table 3. Further, for these solutions are plotted in Fig. 10. We observe that the present method is accurate and even converges for very small values of while FAM and NPCM both diverge. Further, it is noticed that FAM, NPCM, and the present method take execution times of 94.83, 43.56, and 84.57 s respectively. Hence the present method takes less computational time as compared to FAM.

Table 3

Solutions of the system (21) at

Step size	μ	FAM	NPCM	Present method(10)	Exact
10^-2	0.0001	diverges	diverges	4.00001237126	4.0
	0.0005	diverges	diverges	4.00006184560	4.0
10^-3	0.0001	diverges	diverges	4.00001245781	4.0
	0.0005	diverges	diverges	4.00006227827	4.0

Fig. 10

(21)

Example 3

Consider the following fractional order biological-model [38]:

where and . We apply the proposed method to solve the system (22) numerically. For with step-size ; the numerical solutions are plotted in Figs. 3 and 4 individually. Whereas, its phase portraits in versus plane are portrayed in Figs. 5 and 6 separately. Moreover, we found that these graphs match with those obtained by FAM and NPCM reported in [31] and hence validate the applicability of the proposed method.

Fig. 3

(22)

Fig. 4

(22)

Fig. 5

(22)

Fig. 6

(22)

Example 4

Consider the following fractional-order delay system [39]

where and . We perform the numerical simulations for the system (23) using the proposed method. Numerical simulations evidence that this system shows chaotic behavior for and . Further, the numerical solutions (t versus ) and chaotic phase-portraits ( versus ) and ( versus ) are plotted in Figs. 7, 8 and 9 separately. Besides, for and , the stable trajectories and orbits are depicted in Figs. 11 and 12 respectively.

Fig. 7

(23)

Fig. 8

(23)

Fig. 9

(23)

Fig. 11

(23)

Fig. 12

(23)

Example 5

Consider the following fractional-order logistic DDE:

We solve the fractional order system of DDEs (24) using the proposed method. This system shows stable behavior for ()=(0.5,70,26,0.90), periodic oscillations for ()=(0.5,79.3,26,0.90) and chaotic behavior for ()=(0.5, 104, 26, 0.90). In each case, we take step-size . All these solutions are represented graphical in Figs. 13, 14, 15, 16, 17 and 18. Moreover, these graphs match with those obtained in [40].

Fig. 13

(24)

Fig. 14

(24)

Fig. 15

(24)

Fig. 16

(24)

Fig. 17

(24)

Fig. 18

(24)

Conclusions

In this paper, an efficient numerical scheme is developed for solving fractional delay differential equations. Further, the related error analysis of the proposed method is established. Various non-trivial systems of fractional delay differential equations including fractional-order biological models and logistic equations and some chaotic and non-chaotic systems are solved using the proposed method. The absolute and relative errors obtained by FAM, NPCM, and the proposed method are compared numerically and as well as graphically. Numerical simulations show that the proposed method is more accurate than FAM and NPCM and takes less computational time than FAM. Besides, we notice that the proposed method converges even for very small values of the order of fractional derivative operator , when FAM and NPCM diverge.