Numerical Studies for Fractional-Order Logistic Differential Equation with Two Different Delays

Abstract

A numerical method for solving the fractional-order logistic differential equation with two different delays (FOLE) is considered. The fractional derivative is described in the Caputo sense. The proposed method is based upon Chebyshev approximations. The properties of Chebyshev polynomials are utilized to reduce FOLE to a system of algebraic equations. Special attention is given to study the convergence and the error estimate of the presented method. Numerical illustrations are presented to demonstrate utility of the proposed method. Chaotic behavior is observed and the smallest fractional order for the chaotic behavior is obtained. Also, FOLE is studied using variational iteration method (VIM) and the fractional complex transform is introduced to convert fractional Logistic equation to its differential partner, so that its variational iteration algorithm can be simply constructed. Numerical experiment is presented to illustrate the validity and the great potential of both proposed techniques.

1. Introduction

It is known that delay differential equation (DDE) provides a mathematical model for many systems in different fields such as physical, biological systems in which the rate of change of the system depends upon their past history [1, 2]. Introduction of delay in the model enriches its dynamics and allows a precise description of the real life phenomena. DDEs are proved useful in control systems [3], lasers, traffic models [4], metal cutting, epidemiology, neuroscience, population dynamics [2], chemical kinetics [1], and so forth. In DDE one has to provide history of the system over the delay interval [−𝜏;0] as the initial condition. Due to this reason delay systems are infinite dimensional in nature. Because of infinite dimensionality the DDEs are difficult to analyse analytically [5] and hence the numerical solutions play an important role.

On the other side, ordinary and partial fractional differential equations have been the focus of many studies due to their frequent appearance in various applications in fluid mechanics, viscoelasticity, biology, physics and engineering [6]. Consequently, considerable attention has been given to the solutions of fractional differential equations of physical interest. Most fractional differential equations do not have exact solutions, so approximate and numerical techniques [7–12], must be used. Recently, several numerical methods to solve the fractional differential equations have been given such as variational iteration method [13–19], homotopy perturbation method [20], Adomian's decomposition method [21], homotopy analysis method [22], fractional complex transform [15, 23, 24], and collocation method [25–31].

Let us describe some necessary definitions and mathematical preliminaries of the fractional calculus theory required for our subsequent development.

Definition 1.1. The Caputo fractional derivative operator 𝐷^𝛼 of order 𝛼 is defined in the following form: 𝐷^𝛼𝑓(𝑥)=1Γ(𝑚−𝛼)^𝑥₀𝑓^(𝑚)(𝑡)(𝑥−𝑡)^{𝛼−𝑚+1}𝑑𝑡,𝛼>0, (1.1) where 𝑚−1<𝛼≤𝑚,𝑚∈ℕ,𝑥>0.
Similar to integer-order differentiation, Caputo fractional derivative operator is a linear operation: 𝐷^𝛼𝜆𝑓(𝑥)+𝜇𝑔(𝑥)=𝜆𝐷^𝛼𝑓(𝑥)+𝜇𝐷^𝛼𝑔(𝑥), (1.2) where 𝜆 and 𝜇 are constants. For the Caputo's derivative we have [32]: 𝐷^𝛼𝐷𝐶=0,𝐶isaconstant,_𝛼𝑥^𝑛=0,for𝑛∈ℕ₀⌈𝛼⌉;Γ(),𝑛<𝑛+1Γ()𝑥𝑛+1−𝛼^𝑛−𝛼,for𝑛∈ℕ₀,𝑛≥⌈𝛼⌉. (1.3) We use the ceiling function ⌈𝛼⌉ to denote the smallest integer greater than or equal to 𝛼, and ℕ₀={0,1,2,…}. Recall that for 𝛼∈ℕ, the Caputo differential operator coincides with the usual differential operator of integer order.
For more details on fractional derivatives definitions and its properties see [32–34].
In this paper, we consider FOLE with two delays of the form: 𝑑^𝛼𝑥(𝑡)𝑑𝑡^𝛼=𝜌𝑥𝑡−𝑟₁1−𝑥𝑡−𝑟₂,𝑡>0,𝜌>0, (1.4) the parameter 𝛼 refers to the fractional order of time derivative with 0<𝛼≤1.
We also assume an initial condition: 𝑥(𝑡)=𝑥⁰,𝑥⁰>0,𝑡≤0. (1.5) This work is motivated by work done by El-Sayed et al. [8]. The stability of the solution, the existence and the uniqueness of the proposed problem (1.4) are introduced in details in [8].
The main idea of the present work is to apply the Chebyshev collocation method to discretize (1.4) to get a system of algebraic equations thus greatly simplifying the problem. Chebyshev polynomials are a well-known family of orthogonal polynomials on the interval [−1,1] that have many applications [28, 35]. They are widely used because of their good properties in the approximation of functions. However, with our best knowledge, very little work was done to adapt these polynomials to the solution of fractional differential equations.
Khader [25] introduced a new approximate formula of the fractional derivation 𝐷^𝛼𝑥(𝑡) and used it to solve numerically the fractional diffusion equation. In this paper, we will extend this formula to solve the fractional-order logistic equation with two delays (1.4) and prove the error estimate of the introduced formula.
The organization of this paper is as follows. In the next section, the approximation of fractional derivative 𝐷^𝛼𝑥(𝑡) is obtained. Section 3 summarizes the application of Chebyshev collocation method to solve (1.4). In Section 4, the procedure of solution using VIM is given. In Section 5, the use of fractional complex transform is introduced with respect to VIM. Also a conclusion is given in Section 6.

2. Derivation an Approximate Formula for Fractional Derivatives Using Chebyshev Series Expansion

The well-known Chebyshev polynomials [35] are defined on the interval [−1,1] and can be determined with the aid of the following recurrence formula: 𝑇_𝑛+1(𝑧)=2𝑧𝑇_𝑛(𝑧)−𝑇_𝑛−1(𝑧),𝑇₀(𝑧)=1,𝑇₁(𝑧)=𝑧,𝑛=1,2,…. (2.1) The analytic form of the Chebyshev polynomials 𝑇_𝑛(𝑧) of degree 𝑛 is given by the following: 𝑇_𝑛(𝑧)=𝑛^[𝑛/2]_𝑖=0(−1)^𝑖2^{𝑛−2𝑖−1}(𝑛−𝑖−1)!(𝑖)!()!𝑧𝑛−2𝑖^𝑛−2𝑖, (2.2) where [𝑛/2] denotes the integer part of 𝑛/2. The orthogonality condition is ¹₋₁𝑇_𝑖(𝑧)𝑇_𝑗(𝑧)√1−𝑧²𝜋𝑑𝑧=𝜋,for𝑖=𝑗=0;2,for𝑖=𝑗≠0;0,for𝑖≠𝑗. (2.3) In order to use these polynomials on the interval [0,𝐿] we define the so called shifted Chebyshev polynomials by introducing the change of variable 𝑧=(2/𝐿)𝑡−1. The shifted Chebyshev polynomials are defined as 𝑇^∗_𝑛(𝑡)=𝑇_𝑛((2/𝐿)𝑡−1)=𝑇_2𝑛(𝑡/𝐿).

The analytic form of the shifted Chebyshev polynomial 𝑇^∗_𝑛(𝑡) of degree 𝑛 is given by the following: 𝑇^∗_𝑛(𝑡)=𝑛^𝑛_𝑘=0(−1)^𝑛−𝑘2^2𝑘(𝑛+𝑘−1)!𝐿^𝑘()!()!𝑡2𝑘𝑛−𝑘^𝑘. (2.4) The function 𝑥(𝑡), which belongs to the space of square integrable in [0,𝐿], may be expressed in terms of shifted Chebyshev polynomials as follows: 𝑥(𝑡)=^∞_𝑖=0𝑐_𝑖𝑇^∗_𝑖(𝑡), (2.5) where the coefficients 𝑐_𝑖 are given by the following: 𝑐₀=1𝜋^𝐿₀𝑥(𝑡)𝑇^∗₀(𝑡)√𝐿𝑡−𝑡²𝑑𝑡,𝑐_𝑖=2𝜋^𝐿₀𝑥(𝑡)𝑇^∗_𝑖(𝑡)√𝐿𝑡−𝑡²𝑑𝑡,𝑖=1,2,…. (2.6) In practice, only the first (𝑚+1)-terms shifted Chebyshev polynomials are considered. Then we have the following: 𝑥_𝑚(𝑡)=^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖(𝑡). (2.7)

Theorem 2.1 (Chebyshev truncation theorem). The error in approximating 𝑥(𝑡) by the sum of its first 𝑚 terms is bounded by the sum of the absolute values of all the neglected coefficients. If the following: 𝑥_𝑚(𝑡)=^𝑚_𝑘=0𝑐_𝑘𝑇_𝑘(𝑡), (2.8) then 𝐸_𝑇(𝑚)≡𝑥(𝑡)−𝑥_𝑚(𝑡)≤^∞_𝑘=𝑚+1𝑐_𝑘, (2.9) for all 𝑥(𝑡), all 𝑚, and all 𝑡∈[−1,1].

Proof. The Chebyshev polynomials are bounded by one, that is, |𝑇_𝑘(𝑡)|≤1 for all 𝑡∈[−1,1] and for all 𝑘. This implies that the 𝑘-th term is bounded by |𝑐_𝑘|. Subtracting the truncated series from the infinite series, bounding each term in the difference, and summing the bounds gives the theorem.

The main approximate formula of the fractional derivative of 𝑥(𝑡) is given in the following theorem.

Theorem 2.2. Let 𝑥(𝑡) be approximated by Chebyshev polynomials as (2.7) and also suppose 𝛼>0, then: 𝐷^𝛼𝑥_𝑚(𝑡)=^𝑚_⌈𝛼⌉^𝑖_{𝑖=⌈𝛼⌉𝑘=}𝑐_𝑖𝑤^(𝛼)_𝑖,𝑘𝑡^𝑘−𝛼, (2.10) where 𝑤^(𝛼)_𝑖,𝑘 is given by the following: 𝑤^(𝛼)_𝑖,𝑘=(−1)^𝑖−𝑘2^2𝑘𝑖(𝑖+𝑘−1)!Γ(𝑘+1)𝐿^𝑘()!()().𝑖−𝑘2𝑘!Γ𝑘+1−𝛼 (2.11)

Proof. Since the Caputo's fractional differentiation is a linear operation we have 𝐷^𝛼𝑥_𝑚(𝑡)=^𝑚_𝑖=0𝑐_𝑖𝐷^𝛼𝑇^∗_𝑖(𝑡). (2.12) Employing (1.3), in formula (2.4) we have 𝐷^𝛼𝑇^∗_𝑖(𝑡)=0,𝑖=0,1,...,⌈𝛼⌉−1,𝛼>0. (2.13) Also, for 𝑖=⌈𝛼⌉,…,𝑚, and by using (1.3), in formula (2.4) we get 𝐷^𝛼𝑇^∗_𝑖(𝑡)=𝑖^𝑖_{⌈𝛼⌉𝑘=}(−1)^𝑖−𝑘2^2𝑘(𝑖+𝑘−1)!𝐿^𝑘()!()!𝐷𝑖−𝑘2𝑘^𝛼𝑡^𝑘=𝑖^𝑖_{⌈𝛼⌉𝑘=}(−1)^𝑖−𝑘2^2𝑘(𝑖+𝑘−1)!Γ(𝑘+1)𝐿^𝑘()!()()𝑡𝑖−𝑘2𝑘!Γ𝑘+1−𝛼^𝑘−𝛼. (2.14) A combination of (2.13), (2.14), and (2.11) leads to the desired result.

Test example
Consider the function 𝑥(𝑡)=𝑡² with 𝑚=3 and 𝛼=1.5, the shifted Chebyshev series of 𝑡² is 𝑡²=38𝑇^∗₀(𝑡)+48𝑇^∗₁(𝑡)+18𝑇^∗₂(𝑡)+0𝑇^∗₃(𝑡). (2.15) Now, by using (2.10), we obtain 𝐷^3/2𝑡²=^3𝑖_𝑖=2𝑘=2𝑐_𝑖𝑤^3/2_𝑖,𝑘𝑡^𝑘−3/2, (2.16) where, 𝑤^3/2_2,2=16Γ(3/2),𝑤^(3/2)_3,2=−96Γ(3/2),𝑤^(3/2)_3,3=192Γ(),5/2 (2.17) therefore, 𝐷^3/2𝑡²=𝑐₂𝑤^3/2_2,2𝑡^1/2+𝑐₃𝑤^3/2_3,2𝑡^1/2+𝑐₃𝑤^(3/2)_3,3𝑡^3/2=2Γ()𝑡3/2^1/2. (2.18)

Theorem 2.3. The Caputo fractional derivative of order 𝛼 for the shifted Chebyshev polynomials can be expressed in terms of the shifted Chebyshev polynomials themselves in the following form: 𝐷^𝛼𝑇^∗_𝑖(𝑡)=^𝑖_{⌈𝛼⌉𝑘=}^{𝑘−⌈𝛼⌉}_𝑗=0Θ_{𝑖,𝑗,𝑘}𝑇^∗_𝑗(𝑡), (2.19) where Θ_{𝑖,𝑗,𝑘}=(−1)^𝑖−𝑘2𝑖(𝑖+𝑘−1)!Γ(𝑘−𝛼+1/2)𝐿^𝑘−𝛼ℎ_𝑗Γ(𝑘+1/2)(𝑖−𝑘)!Γ𝑘−𝛼−𝑗+1Γ𝑘+𝑗−𝛼+1,𝑗=0,1,…. (2.20)

Proof. We concern the properties of the shifted Chebyshev polynomials [35] and expand 𝑡^𝑘−𝛼in (2.14) in the following form [36]: 𝑡^𝑘−𝛼=^{𝑘−⌈𝛼⌉}_𝑗=0𝑐_𝑘𝑗𝑇^∗_𝑗(𝑡), (2.21) where 𝑐_𝑘𝑗 can be obtained using (2.6), where 𝑥(𝑡)=𝑡^𝑘−𝛼 then, 𝑐_𝑘𝑗=2ℎ_𝑗𝜋^𝐿₀𝑡^𝑘−𝛼𝑇^∗_𝑗(𝑡)√𝐿𝑡−𝑡²𝑑𝑡,ℎ₀=2,ℎ_𝑗=1,𝑗=1,2,…. (2.22) At 𝑗=0 we find, 𝑐_𝑘0=(1/𝜋)∫^𝐿₀(𝑡^𝑘−𝛼𝑇^∗₀(𝑡)/√𝐿𝑡−𝑡²)𝑑𝑡=(𝐿^𝑘−𝛼/√𝜋)(Γ(𝑘−𝛼+1/2)/Γ(𝑘−𝛼+1)), also, at any 𝑗 and using the formulae (2.4) and (2.6) we can find that 𝑐_𝑘𝑗=𝑗√𝜋^𝑗_𝑟=0(−1)^𝑗−𝑟𝑗+𝑟−1!2^2𝑟+1Γ(𝑘+𝑟−𝛼+1/2)𝐿^𝑘−𝛼𝑗−𝑟!(2𝑟)!Γ(𝑘+𝑟−𝛼+1),𝑗=1,2,3,…, (2.23) employing (2.14) and (2.21) gives 𝐷^𝛼𝑇^∗_𝑖(𝑡)=^𝑖_{⌈𝛼⌉𝑘=}^{𝑘−⌈𝛼⌉}_𝑗=0Θ_{𝑖,𝑗,𝑘}𝑇^∗_𝑗(𝑡),𝑖=⌈𝛼⌉,⌈𝛼⌉+1,…, (2.24) where Θ_{𝑖,𝑗,𝑘}=𝑖(−1)^𝑖−𝑘(𝑖+𝑘−1)!2^2𝑘𝑘!Γ(𝑘−𝛼+1/2)𝐿^𝑘−𝛼(𝑖−𝑘)!(2𝑘)!√𝜋(Γ(𝑘+1−𝛼))²()𝑗=0;−1^𝑖−𝑘𝑖𝑗(𝑖+𝑘−1)!2^2𝑘+1𝑘!√()!()!×𝜋Γ𝑘+1−𝛼)(𝑖−𝑘2𝑘_𝑗𝑟∑₌₀(−1)^𝑗−𝑟𝑗+𝑟−1!2^2𝑟Γ(𝑘+𝑟−𝛼+1/2)𝐿^𝑘−𝛼𝑗−𝑟!(2𝑟)!Γ(𝑘+𝑟−𝛼+1),𝑗=1,2,3,…. (2.25) After some lengthly manipulation Θ_{𝑖,𝑗,𝑘} can put in the following form: Θ_{𝑖,𝑗,𝑘}=(−1)^𝑖−𝑘2𝑖(𝑖+𝑘−1)!Γ(𝑘−𝛼+1/2)𝐿^𝑘−𝛼ℎ_𝑗Γ(𝑘+1/2)(𝑖−𝑘)!Γ𝑘−𝛼−𝑗+1Γ𝑘+𝑗−𝛼+1,𝑗=0,1,..., (2.26) and this completes the proof of the theorem.

Theorem 2.4. The error |𝐸_𝑇(𝑚)|=|𝐷^𝛼𝑥(𝑡)−𝐷^𝛼𝑥_𝑚(𝑡)| in approximating 𝐷^𝛼𝑥(𝑡) by 𝐷^𝛼𝑥_𝑚(𝑡) is bounded by the following: 𝐸_𝑇(𝑚)≤^∞_𝑖=𝑚+1𝑐_𝑖^𝑖_{⌈𝛼⌉𝑘=}^{𝑘−⌈𝛼⌉}_𝑗=0Θ_{𝑖,𝑗,𝑘}. (2.27)

Proof. A combination of (2.5), (2.7), and (2.19) leads to the following: 𝐸_𝑇(𝑚)=𝐷^𝛼𝑥(𝑡)−𝐷^𝛼𝑥_𝑚(𝑡)=^∞_𝑖=𝑚+1𝑐_𝑖^𝑖_{⌈𝛼⌉𝑘=}^{𝑘−⌈𝛼⌉}_𝑗=0Θ_{𝑖,𝑗,𝑘}𝑇^∗_𝑗(𝑡), (2.28) but |𝑇^∗_𝑗(𝑡)|≤1, so, we can obtain 𝐸_𝑇(𝑚)≤^∞_𝑖=𝑚+1𝑐_𝑖^𝑖_{⌈𝛼⌉𝑘=}^{𝑘−⌈𝛼⌉}_𝑗=0Θ_{𝑖,𝑗,𝑘}, (2.29) and subtracting the truncated series from the infinite series, bounding each term in the difference, and summing the bounds completes the proof of the theorem.

3. Procedure of Solution for the Fractional-Order Logistic Equation

Consider the FOLE of type given in (1.4). In order to use the Chebyshev collocation method, we first approximate 𝑥(𝑡) as follows: 𝑥_𝑚(𝑡)=^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖(𝑡). (3.1) From (1.4), (3.1), and Theorem 2.2 we have ^𝑚_⌈𝛼⌉^𝑖_{𝑖=⌈𝛼⌉𝑘=}𝑐_𝑖𝑤^(𝛼)_𝑖,𝑘𝑡^𝑘−𝛼=𝜌^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖𝑡−𝑟₁1−^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖𝑡−𝑟₂. (3.2) We now collocate (3.2) at (𝑚+1−⌈𝛼⌉) points 𝑡_𝑝 as follows: ^𝑚_⌈𝛼⌉^𝑖_{𝑖=⌈𝛼⌉𝑘=}𝑐_𝑖𝑤^(𝛼)_𝑖,𝑘𝑡_𝑝^𝑘−𝛼=𝜌^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖𝑡_𝑝−𝑟₁1−^𝑚_𝑖=0𝑐_𝑖𝑇^∗_𝑖𝑡_𝑝−𝑟₂,𝑝=0,1,...,𝑚−⌈𝛼⌉. (3.3) For suitable collocation points we use roots of shifted Chebyshev polynomial 𝑇^∗_{𝑚+1−⌈𝛼⌉}(𝑡). Also, by substituting (2.7) in the initial condition (1.5) we can find the following equation: ^𝑚_𝑖=0(−1)^𝑖𝑐_𝑖=𝑥⁰. (3.4) Equation (3.3), together with the equation of the initial condition (3.4), give (𝑚+1) of nonlinear algebraic equations which can be solved, for the unknown 𝑐_𝑖,𝑖=0,1,…,𝑚.

In the following section, to achieve from the validity and the accuracy we compare our approximate solution with those obtained using the variational iteration method.

4. Procedure of Solution Using VIM

The VIM gives the possibility to write the solution of (1.4), (0<𝛼≤1) with the aid of the correction functionals in the form: 𝑥_𝑛+1(𝑡)=𝑥_𝑛(𝑡)+^𝑡₀𝜆(𝜏)𝑑𝑥_𝑛𝑥𝑑𝜏−𝜌_𝑛𝑡−𝑟₁1−𝑥_𝑛𝑡−𝑟₂𝑑𝜏, (4.1) where 𝜆 is the general Lagrange multiplier, which can be identified optimally via the variational theory [13]. The function 𝑥_𝑛 is a restricted variation, which means that 𝛿𝑥_𝑛=0. Therefore, we first determine the Lagrange multiplier 𝜆 that will be identified optimally via integration by parts. The successive approximations 𝑥_𝑛,𝑛≥0, of the solution 𝑥(𝑡) will be readily obtained upon using the obtained Lagrange multiplier and by using any selective function 𝑥₀. The initial values of the solution are usually used for selecting the zeroth approximation 𝑥₀. With 𝜆 determined, several approximations 𝑥_𝑛,𝑛≥0, follow immediately. Making the above correction functional stationary: 𝛿𝑥_𝑛+1(𝑡)=𝛿𝑥_𝑛(𝑡)+𝛿^𝑡₀𝜆(𝜏)𝑑𝑥_𝑛𝑥𝑑𝜏−𝜌_𝑛𝑡−𝑟₁1−𝑥_𝑛𝑡−𝑟₂𝑑𝜏=𝛿𝑥_𝑛(𝑡)+^𝑡₀𝛿𝜆(𝜏)𝑑𝑥_𝑛𝑑𝜏𝑑𝜏=𝛿𝑥_𝑛(𝑡)+𝜆(𝜏)𝛿𝑥_𝑛(𝜏)_𝜏=𝑡−^𝑡₀𝛿𝑥_𝑛𝜆̇(𝜏)𝑑𝜏=0, (4.2) where 𝛿𝑥_𝑛 is considered as a restricted variation, that is, 𝛿𝑥_𝑛=0, yields the following stationary conditions (by comparison the two sides in the above equation): 𝜆̇(𝜏)=0,1+𝜆(𝜏)|_𝜏=𝑡=0. (4.3) The equation in (4.3) is called Lagrange-Euler equation with the natural boundary condition. The solution of this equation gives the Lagrange multiplier 𝜆(𝜏)=−1.

Now, by substituting in (4.1), we obtained 𝑥_𝑛+1(𝑡)=𝑥_𝑛(𝑡)−^𝑡₀𝑑^𝛼𝑥_𝑛𝑑𝜏^𝛼−𝜌𝑥_𝑛𝑡−𝑟₁1−𝑥_𝑛𝑡−𝑟₂𝑑𝜏,𝑛≥0. (4.4)

5. Procedure of Solution Using VIM with Fractional Complex Transform

In this section, we use the fractional complex transform [23, 24] to convert the fractional Logistic equation into its differential partner, so that the VIM can be effectively used. To achieve this aim we can use the following fractional complex transform: 𝑡𝑇=^𝛼Γ(),1+𝛼 (5.1) and using Jumarie's chain rule [37, 38], we have 𝑑^𝛼𝑥𝑑𝑡^𝛼=𝑑𝑥⋅𝑑𝑑𝑇^𝛼𝑇𝑑𝑡^𝛼=𝑑𝑥.𝑑𝑇 (5.2) So, we can obtain the corresponding ordinary differential equation of (1.4) as follows: 𝑑𝑥(𝑇)𝑑𝑇−𝜌𝑥𝑇−𝑟₁1−𝑥𝑇−𝑟₂=0,𝑇>0,𝜌>0. (5.3) The variational iteration algorithm is 𝑥_𝑛+1(𝑇)=𝑥_𝑛(𝑇)−^𝑇₀𝑑𝑥(𝑠)𝑑𝑠−𝜌𝑥𝑠−𝑟₁1−𝑥𝑠−𝑟₂𝑑𝑠. (5.4) We start with initial approximation 𝑥₀(𝑡)=0.5, and by using the above iteration formula (4.4), we can directly obtain the components of the solution. Consequently, the exact solution may be obtained by using the following: 𝑥(𝑡)=lim_𝑛→∞𝑥_𝑛(𝑡). (5.5) Now, the first two components of the solution 𝑥(𝑡) at 𝛼=1,𝜌=0.5 of the fractional Logistic equation with different two delays by using (4.4) or (5.4) are 𝑥₀(𝑡)=0.5,  𝑥₁(𝑡)=0.5+𝑡⁴+0.0555556𝑡⁹,…. The numerical results of this example are presented in Figures 1–5, using the Chebyshev collocation method (ChebM) and compared with the approximate solution using VIM. We considered different values of the delay parameters 𝑟₁ and 𝑟₂ with different values of 𝑚 and 𝛼.

(a)

(b)

Figure 1 

The behavior of the numerical solution (at 𝛼=1,𝜌=0.5) using ChebM with (𝑚=6) and VIM at 𝑟1=0.0, 𝑟2=0.7 (a) and 𝑟1=𝑟2=0.7 (b).

Figure 2 

The behavior of the numerical solution (at 𝛼=1,𝜌=0.5) using ChebM with (𝑚=6) and VIM at 𝑟1=0.7, 𝑟2=0.7.

(a)

(b)

Figure 3 

The behavior of the numerical solution (at 𝛼=0.85,𝜌=0.5) using ChebM with (𝑚=6) at 𝑟1=0.0, 𝑟2=0.7 (a) and 𝑟1=𝑟2=0.7 (b).

(a)

(b)

Figure 4 

Show phase portrait of the system, that is, plot of 𝑥(𝑡) versus 𝑥(𝑡−𝑟1) (a) and plot of 𝑥(𝑡) versus 𝑥(𝑡−𝑟2) (b) at 𝑟1=0.7, 𝑟2=0.7, and 𝛼=0.85.

(a)

(b)

Figure 5 

Show phase portrait of the system, that is, plot of 𝑥(𝑡) versus 𝑥(𝑡−𝑟1) (a) and plot of 𝑥(𝑡−𝑟1) versus 𝑥(𝑡−𝑟2) (b) at 𝑟1=0.7, 𝑟2=0.7, and 𝛼=0.85.

In Figure 1, the behavior of the numerical solution (at 𝛼=1,𝜌=0.5) using ChebM with (𝑚=6) and VIM at 𝑟₁=0.0, 𝑟₂=0.7 (a) and 𝑟₁=𝑟₂=0.7 (b) is presented in the interval [0,1]. In Figure 2, the behavior of the numerical solution (at 𝛼=1,𝜌=0.5) using ChebM with (𝑚=6) and VIM at 𝑟₁=0.7, 𝑟₂=0.7 is presented in the interval [0,2]. In Figure 3, the behavior of the numerical solution (at 𝛼=0.85,𝜌=0.5) using ChebM with (𝑚=6) at 𝑟₁=0.0, 𝑟₂=0.7 (a) and 𝑟₁=𝑟₂=0.7 (b) is presented in the interval [0,10]. Also, from Figures 4 and 5 we can see that the system shows a periodic (chaotic) behavior. In the above experiments we have decreased the value of 𝛼 and observed that the system becomes periodic for 𝛼.

In Table 1, we presented the CPU time needed for the computation with each method using different values of 𝑚, where 𝑚 represents the iteration number of VIM and represents the number of term of the series using the ChebM.

Table 1 

The CPU time needed for the computation with each method.

From Figure 1, we can conclude that the obtained numerical results using the proposed method are in excellent agreement with the exact solution and the solution using VIM. On the other hand, in Figure 2, we saw that the solution using VIM is not agreement with the exact solution at a large domain. This result confirm that our proposed method is more accurate and efficiency from VIM.

6. Conclusion

In this paper, some interesting fractional delay differential equations arising in biology have been solved. It is observed that even two dimensional delayed systems of fractional order show chaotic behavior, and below some critical order, the system changes its nature and becomes periodic. In some cases it is observed that the phase portrait gets stretched as the order of the derivative is reduced. We used two computational methods, Chebyshev spectral method and variational iteration method for solving the fractional-order Logistic equation with two different delays. We derived an approximate formula of the fractional derivative. The properties of the Chebyshev polynomials are used to reduce FOLE with two different delays to the solution of nonlinear system of algebraic equations. Special attention is given to study the convergence analysis and estimate the upper bound of the error of the derived formula. From the solutions obtained using the suggested method we can conclude that these solutions are in excellent agreement with the exact solution and show that these approaches can solve the problem effectively. It is evident that the overall errors can be made smaller by adding new terms from the series (2.7). Comparisons are made between approximate solutions to illustrate the validity and the great potential of the proposed techniques. All numerical results are obtained using MatLab 8.