# Numerical Studies for Fractional Functional Differential Equations with Delay Based on BDF-Type Shifted Chebyshev Approximations

Authors: V. G. Pimenov; A. S. Hendy
Journal: Abstract and Applied Analysis
 (2015)

---

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


---

## Body

## 1. Introduction

Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences [1–4]. Some approximate and numerical methods are used to obtain approximate solutions [5–11]. The general theory of differential equations with delays (DDEs) is widely developed and discussed in the literature [12–16]. DDEs are differential equations in which the rate of change of y(t) does not depend only on the values of y for the same time value but also on time values less than t. In the simplest case, DDEs have the form(1)y′t=ft,yt,yt-τ,t∈0,L,τ>0,under the following initial conditions: (2)y(t)=ϕ(t),-τ≤t≤0,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 [19–23]. Some approaches are built to obtain numerical and approximate solutions for FDDEs:(3)Dβyt=ft,yt,yt-τ,t∈0,L,t≥0,m-1<β≤m,m=1,2,3,…,under the following initial conditions: (4)y(t)=ϕ(t),t≤0.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 [25–27]. Bhalekar and Daftardar-Gejji adapted the Adams Bashforth Moulton predictor corrector method (EABMPC method) [28] to solve this form of FDDEs (3)-(4), where 0<β≤1. 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: (5)Dβyt=ayt-τ+byt+ft,t≥0,0<β≤1,under the following initial conditions: (6)y(t)=ϕ(t),t∈[-τ,0].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:(7)Dβyt=ft,yt,yt-τ,t∈0,L,t≥0,m-1<β≤m,y(t)=ϕ(t),t∈[-τ,0].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 f 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 Section2, 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 nβ such that n 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 operatorD(ν) of order ν is defined in the following form: (8)D(ν)f(x)=1Γ(m-ν)∫0x‍f(m)(ξ)x-ξν-m+1dξ,ν>0,where m-1<ν≤m, m∈N, x>0.Similar to integer-order differentiation, the Caputo fractional derivative operator is a linear operation:(9)D(ν)(λp(x)+μq(x))=λD(ν)p(x)+μD(ν)q(x),where λ and μ are constants.For Caputo’s derivative we have(10)D(ν)C=0,Cisaconstant,D(ν)xn=0,forn∈N0,n<ν;Γ(n+1)Γ(n+1-ν)xn-ν,forn∈N0,n≥ν.We use the ceiling function ν to denote the smallest integer greater than or equal to ν and N0={0,1,2,…}. Recall that, for ν∈N, the Caputo differential operator coincides with the usual differential operator of integer order.For more details on fractional derivatives definitions and their properties, see [35–38].


### 2.2. The Definition and Properties of the Shifted Chebyshev Polynomials

The well-known Chebyshev polynomials [39] are defined on the interval [-1,1] and can be determined with the aid of the following recurrence formula:(11)Tn+1z=2zTnz-Tn-1z,T0z=1,T1z=z,n=1,2,….It is well known that Tn(1)=1, Tn(-1)=(-1)n. The analytic form of the Chebyshev polynomials Tn(z) of degree n is given by (12)Tn(z)=∑i=0n/2‍-1i2n-2i-1n(n-i-1)!(i)!(n-2i)!zn-2i,where n/2 denotes the integer part of n/2. The orthogonality condition is(13)∫-11‍Ti(z)Tj(z)1-z2dz=π,fori=j=0;π2,fori=j≠0;0,fori≠j.In order to use these polynomials on the interval [0,L] we define the so-called shifted Chebyshev polynomials by introducing the change of variable z=2x/L-1. The shifted Chebyshev polynomials are defined as(14)Tn∗(x)=Tn(2xL-1),whereT0∗(x)=1,T1∗(x)=2xL-1.The analytic form of the shifted Chebyshev polynomial Tn∗(x) of degree n is given by (15)Tn∗(x)=n∑k=0n‍-1n-k(n+k-1)!22k(n-k)!(2k)!Lkxk,where Tn∗(0)=(-1)n and Tn∗(L)=1. The orthogonality condition of these polynomials is (16)∫0L‍Tj∗(x)Tk∗(x)w(x)dx=δjkhk,where the weight function w(x)=1/Lx-x2, hk=(bk/2)π, with b0=2, bk=1, k≥1.The functiony(x), which belongs to the space of square integrable in [0,L], may be expressed in terms of shifted Chebyshev polynomials as (17)y(x)=∑n=0∞‍cnTn∗(x),where the coefficients cn are given by (18)cn=1hn∫0L‍yxTn∗xwxdx,n=0,1,2,….The well-known shifted Chebyshev polynomials of the first kind of degree n are defined on the interval [0,L] as in (15). We choose the grid (interpolation) points to be the Chebyshev-Gauss Lobatto points associated with the interval [0,L], xr=L/2-(L/2)cos⁡(πr/N), r=0,1,…,N. These grids can be written as L=xN<xN-1<⋯<x1<x0=0.Clenshaw and Curtis [40] introduced an approximation of the function y(x), and we reformulate it to be used on the shifted Chebyshev polynomials as follows: (19)yN(x)=∑n=0N′′anTn∗(x),an=2N∑r=0N‍′′y(xr)Tn∗(xr).The summation symbol with double primes denotes a sum with both first and last terms halved.


## 2.1. The Caputo Fractional Derivatives

Definition 1.
The Caputo fractional derivative operatorD(ν) of order ν is defined in the following form: (8)D(ν)f(x)=1Γ(m-ν)∫0x‍f(m)(ξ)x-ξν-m+1dξ,ν>0,where m-1<ν≤m, m∈N, x>0.Similar to integer-order differentiation, the Caputo fractional derivative operator is a linear operation:(9)D(ν)(λp(x)+μq(x))=λD(ν)p(x)+μD(ν)q(x),where λ and μ are constants.For Caputo’s derivative we have(10)D(ν)C=0,Cisaconstant,D(ν)xn=0,forn∈N0,n<ν;Γ(n+1)Γ(n+1-ν)xn-ν,forn∈N0,n≥ν.We use the ceiling function ν to denote the smallest integer greater than or equal to ν and N0={0,1,2,…}. Recall that, for ν∈N, the Caputo differential operator coincides with the usual differential operator of integer order.For more details on fractional derivatives definitions and their properties, see [35–38].


## 2.2. The Definition and Properties of the Shifted Chebyshev Polynomials

The well-known Chebyshev polynomials [39] are defined on the interval [-1,1] and can be determined with the aid of the following recurrence formula:(11)Tn+1z=2zTnz-Tn-1z,T0z=1,T1z=z,n=1,2,….It is well known that Tn(1)=1, Tn(-1)=(-1)n. The analytic form of the Chebyshev polynomials Tn(z) of degree n is given by (12)Tn(z)=∑i=0n/2‍-1i2n-2i-1n(n-i-1)!(i)!(n-2i)!zn-2i,where n/2 denotes the integer part of n/2. The orthogonality condition is(13)∫-11‍Ti(z)Tj(z)1-z2dz=π,fori=j=0;π2,fori=j≠0;0,fori≠j.In order to use these polynomials on the interval [0,L] we define the so-called shifted Chebyshev polynomials by introducing the change of variable z=2x/L-1. The shifted Chebyshev polynomials are defined as(14)Tn∗(x)=Tn(2xL-1),whereT0∗(x)=1,T1∗(x)=2xL-1.The analytic form of the shifted Chebyshev polynomial Tn∗(x) of degree n is given by (15)Tn∗(x)=n∑k=0n‍-1n-k(n+k-1)!22k(n-k)!(2k)!Lkxk,where Tn∗(0)=(-1)n and Tn∗(L)=1. The orthogonality condition of these polynomials is (16)∫0L‍Tj∗(x)Tk∗(x)w(x)dx=δjkhk,where the weight function w(x)=1/Lx-x2, hk=(bk/2)π, with b0=2, bk=1, k≥1.The functiony(x), which belongs to the space of square integrable in [0,L], may be expressed in terms of shifted Chebyshev polynomials as (17)y(x)=∑n=0∞‍cnTn∗(x),where the coefficients cn are given by (18)cn=1hn∫0L‍yxTn∗xwxdx,n=0,1,2,….The well-known shifted Chebyshev polynomials of the first kind of degree n are defined on the interval [0,L] as in (15). We choose the grid (interpolation) points to be the Chebyshev-Gauss Lobatto points associated with the interval [0,L], xr=L/2-(L/2)cos⁡(πr/N), r=0,1,…,N. These grids can be written as L=xN<xN-1<⋯<x1<x0=0.Clenshaw and Curtis [40] introduced an approximation of the function y(x), and we reformulate it to be used on the shifted Chebyshev polynomials as follows: (19)yN(x)=∑n=0N′′anTn∗(x),an=2N∑r=0N‍′′y(xr)Tn∗(xr).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 form(20)Dβyt=ft,yt,yt-τ,t∈0,L,t≥0,yt=ϕt,t∈-τ,0,such that m-1<β≤m and m is an integer.Yang and Cao studied the existence and uniqueness of these types of (20) using the fixed point theory [24]. Now, let ys be an approximation of the theoretical solution y(t) at ts. We are interested in obtaining a numerical approximation at the point ts+1=ts+h such that h is the step size.If we rewrite the solutiony(t) in the interval [ts,ts+1] in terms of a new variable α defined by(21)t=ts+hαL,α∈0,L,we can approximate the solution of FDDEs expressed in the form (22)y(t)=y(ts+hαL)=y-(α).Approximate (22) in terms of Clenshaw and Curtis finite sum approximations as follows: (23)y-(α)=∑n=0N‍′′anTn∗(α),an=2N∑r=0N‍′′y(ts+ζrh)Tn∗(αr),such that (24)ζr=αrL,αr=L2-L2cos⁡πrN,Dt(β)y-(α)=2N∑n=0N‍′′∑r=0N‍′′yts+ζrhTn∗αrDtβTn∗αt.Using the analytic form of shifted Chebyshev polynomials (15) and the properties of Caputo fractional derivatives (10) gives(25)Dt(β)Tn∗(α(t))=∑k=0n‍n-1n-k22kn+k-1!n-k!2k!LkDtβLht-Lhtsk(26)=∑k=βn‍n-1n-k22kn+k-1!Γk+1n-k!2k!LkΓk+1-β·LhβLht-tsk-β(27)=∑k=βn‍n-1(n-k)2(2k)(n+k-1)!Γ(k+1)Lβ(n-k)!(2k)!LkΓ(k+1-β)hβαk-β.Now, αk-β can be expressed approximately in terms of shifted Chebyshev series as follows: (28)αk-β≅∑j=0N‍ckjTj∗(α),where ckj is obtained from (18) as clarified by Doha et al. [5],(29)ckj=1πLk-βΓ(k-β+1/2)Γ(k-β+1),forj=0;jLk-βπ·∑r=0j‍-1j-r(j+r-1)!22r+1Γ(k+r-β+1/2)(j-r)!(2r)!Γ(k+r-β+1),forj=1,2,….From (27)–(29), we get(30)DtβTn∗αt=1hβ∑j=0∞‍ωβn,jTj∗α,n=β,β+1,…,where ωβ(n,j)=∑k=βn‍ψn,j,k,(31)ψn,j,k=n-1n-k(n+k-1)!22kk!Γ(k-β+1/2)Lβ(n-k)!(2k)!πΓk-β+12,forj=0;-1n-knj(n+k-1)!22k+1k!Lβ(n-k)!(2k)!Γ(k-β+1)π·∑r=0j‍-1j-r(j+k-1)!22rΓ(k+r-β+1/2)(j-r)!(2r)!Γ(k+r-β+1),forj=1,2,….Then, ψn,j,k can be put in the following form:(32)ψn,j,k=-1n-k2n(n+k-1)!Γ(k-β+12)·bjLβΓk+12n-k!Γk-β-j+1·Γk+j-β+1Lβ-1,j=0,1,….Then,(33)Dt(β)Tn∗(α(t))=1hβ∑j=0N‍∑k=βn‍-1n-k2n(n+k-1)!Γ(k-β+12)Γk-β-j+1bj-1·bjΓ(k+12)(n-k)!Γ(k-β+j+1)·Γk-β-j+1bj-1Tj∗α.From (23) and (33), we obtain the following approximate equation: (34)f(t,y(t),y(t-τ))=4Nhβ∑n=βN′′∑r=0N‍′′∑j=0N‍∑k=βn‍-1n-kn(n+k-1)!Γ(k-β+1/2)y(ts+ζrh)Tn∗(αr)bjΓ(k+1/2)(n-k)!Γ(k-β+j+1)Γ(k-β-j+1)Tj∗(α).Evaluating the formula (34) at αl=L/2-(L/2)cos⁡(πl/N), l=1,2,…,N gives an implicit system of algebraic equations(35)f(ts+ζlh,y(ts+ζlh),y(ts+ζlh-τ))=4Nhβ∑n=βN′′∑r=0N‍′′∑j=0N‍∑k=βn‍-1n-kn(n+k-1)!Γ(k-β+1/2)y(ts+ζrh)Tn∗(αr)bjΓ(k+1/2)(n-k)!Γ(k-β+j+1)Γ(k-β-j+1)Tj∗(αl).It can be observed that, for any delayed termτ, (ts+ζlh-τ) may not be a grid point tl for any l. So, we establish the approximation for the delayed function y(t-τ) as follows:(36)dl=y(tl-τ(tl))=y(ts+ζlh-τ(ts+ζlh)).So, we define the delayed term as τ=(ms+δl)h such that 0≤δl<1 and ms is a positive integer. As ts+1=ts+h, then (37)dl=y(t0+(s-ms)h+(ζl-δl)h).Then, we have two different cases.The First Case (
δ
l
=
0
). Consider(38)dl=y(ts-ms+ζlh),fors>ms;Λl,fors≤ms,such that (39)Λl=y(t0+ζlh),fors=ms;y(t0),fors+ζl=ms;ϕ(tl-τ(tl)),fors+ζl<ms.The Second Case (
0
<
δ
l
<
1
). In this case, y(t-τ) is not a grid point. To overcome this, we introduce the following approximation:(40)dl=y(t0+(s-ms-1)h+ζlh-(1-δl)h)=y((ts-ms-1+ζlh)+(1-δl)h)≃y(ts-ms-1+ζlh)+1-δlhy′(ts-ms-1+ζlh)+O(h2)≃yts-ms-1+ζlh+1-δlh·yts-ms-1+ζlh+h-yts-ms-1+ζlhh≃(1-δl)y(ts-ms+ζlh)+δly(ts-ms-1+ζlh).Then, we have (41)dl=(1-δl)y(ts-ms+ζlh)+δly(ts-ms-1+ζlh),fors>ms;(1-δl)y(t0+ζlh)+δlϕ(tl-τ(tl)),fors=ms;ħl,fors<ms,such that (42)ħl=(1-δl)y(t0)+δlϕ(tl-τ(tl)),fors+ζl=ms;ϕ(tl-τ(tl)),fors+ζl<ms.From (35)–(42), we have the following numerical approximation scheme:(43)hβf(ts+ζlh,y(ts+ζlh),dl)=4N∑n=βN′′∑r=0N‍′′∑j=0N∑k=βn‍-1n-kn(n+k-1)!Γ(k-β+1/2)y(ts+ζrh)Tn∗(αr)bjΓ(k+1/2)(n-k)!Γ(k-β+j+1)Γ(k-β-j+1)Tj∗(αl),where the unknowns are the values of the solution at the intermediate points y(ts+ζlh).For DDEs(β=1), explicitly, if we put N=4 in (43), the algebraic system will have the following form:(44)hf1=-(2+2)y(ts)+2y(ts+ζ1h)+22y(ts+ζ2h)-2y(ts+ζ3h)+(2-2)y(ts+ζ4h),hf2=1y(ts)-22y(ts+ζ1h)+0y(ts+ζ2h)+22y(ts+ζ3h)-1y(ts+ζ4h),hf3=-(2-2)y(ts)+2y(ts+ζ1h)-22y(ts+ζ2h)-2y(ts+ζ3h)+(2+2)y(ts+ζ4h),hf4=1y(ts)-(8-42)y(ts+ζ1h)+4y(ts+ζ2h)-(8+42)y(ts+ζ3h)+11y(ts+ζ4h).For FDDEs (β=1/2), also, if we put N=4 in (43), the algebraic system will have the following form:(45)h0.5f1=-(2(335+1822)a)y(ts)+280+3202ayts+ζ1h-1221-72ayts+ζ2h-((56+162)a)y(ts+ζ3h)+((-28+432)a)y(ts+ζ4h),h0.5f2=-(712d)y(ts)-((168+402)d)y(ts+ζ1h)+(1802d)y(ts+ζ2h)+((168-402)d)y(ts+ζ3h)-((292)d)8y(ts+ζ4h),h0.5f3=-(2(335-1822)b)y(ts)+((56-162)b)y(ts+ζ1h)-1221+72byts+ζ2h+-280+3202byts+ζ3h+((28+432)b)y(ts+ζ4h),h0.5f4=-(432c)y(ts)+((168-1522)c)y(ts+ζ1h)+(122c)y(ts+ζ2h)-((168+1522)c)y(ts+ζ3h)+3352cyts+ζ4h,such that(46)a=11052πsin(π8),b=11052πcos⁡(π8),c=11052π,d=105π.For the sake of simplicity, we take this abbreviation (47)fl=f(ts+ζlh,y(ts+ζlh),dl),l=1,2,3,4.Solving this system, we obtain in particular the required value at the final point on the interval y(ts+ζ4h)=y(ts+h)=ys+1. If we repeat this procedure along the integration interval [0,L], 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 N=4 can be expressed by one step recurrence formula: (48)CYs=BYs-1+hβFs,where(49)Fs=f1,f2,f3,f4T,Ys=yts+ζ1h,yts+ζ2h,yts+ζ3h,yts+ζ4hT,Ys-1=yts-1+ζ1h,yts-1+ζ2h,yts-1+ζ3h,yts-1+ζ4hT.Remark 2.
Consider(50)y(ts-ζ4-lh)=y(ts-1+ζlh),l=1,2,3,4,Ys-1=yts+ζ3h,yts+ζ2h,yts+ζ1h,yts+ζ0hT.For DDEs (β=1), the constant matrices B and C are given by (51)B=0002+200010002-2000-1,C=222-22-2-22022-12-22-22+2-8+424-8-4211,which have a completely agreement with the work done by Ramos and Vigo-Aguiar [34].For FDDEs(β=0.5), the constant matrices B and C are also given by(52)B=0002(335+1822)a000712d0002(335-1822)b000432c,C=(280+3202)a-122(1-72)a-(56+162)a(-28+432)a(-168-402)d1802d(168-402)d(-292)d(56-162)b-122(1+72)b(-280+3202)b(28+432)b(168-1522)c122c(-168-1522)c3352c.As the matrix C is nonsingular, multiplying the formula (48) by the inverse matrix C-1 yields (53)Ys=C-1BYs-1+hβC-1Fs.Remark 3.
The productC-1B is a square matrix of dimension 4 on the following form:(54)0001000100010001.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(β=1), 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(β=0.5), see Table 2.The contents of Butcher tableau at β=0.5 can be written also in terms of fractions but it will be complicated. So, we write it in decimal form.
In the names of Remark4, we can say that the theory of Runge-Kutta methods may be applied here. For our approach, when we deal with DDEs (β=1), the Butcher conditions for order four are verified. So, this method has the fourth order for DDEs.Table 1
α
1
1
/
384
(
88
-
4
2
)
1
/
384
(
40
-
64
2
)
1
/
384
(
88
-
28
2
)
-
1
/
16
α
2
1
/
384
(
64
+
48
2
)
1
/
6
1
/
384
(
64
-
48
2
)
0
α
3
1
/
384
(
88
+
28
2
)
1
/
384
(
40
+
64
2
)
1
/
384
(
88
+
4
2
)
-
1
/
16
α
4
1
/
3
1
/
3
1
/
3
0
1
/
3
1
/
3
1
/
3
0Table 2
α
1
0.414917
-
0.0779
0.0557512
-
0.0254207
α
2
0.317636
0.528742
-
0.101593
0.0415419
α
3
0.172975
0.420603
0.495111
-
0.0645867
α
4
0.169572
0.313329
0.409383
0.221557
0.169572
0.313329
0.409383
0.221557Consider the difference operators associated with each stage in our method as follows:(55)Liyt,h=yts+ζih-yts-hβ∑j=1N‍cijyβts+ζih,i=1,2,3,4,(56)Li(y(t),h)=y(ts+ζih)-y(ts)-hβ∑j=1N‍cijf(ts+ζih,y(ts+ζih),di),and using the fractional Taylor series about ts, we obtain(57)yts+ζih=y(ts)+yβ(ts)ζihβΓβ+1+y2β(ts)ζih2βΓ2β+1+⋯+y((N+1)β)(ts+θζih)ζih(N+1)βΓ((N+1)β+1),f(ts+ζih,y(ts+ζih),di)=fi(ζih)=fi(0)+fiβ(0)ζihβΓβ+1+fi2β(0)ζih2βΓ2β+1+⋯+fiNβθζihζihNβΓNβ+1.Theorem 5.
If the proposed Runge-Kutta method with shifted Chebyshev polynomials (48) is of order Nβ and if all partial fractional derivatives of f(t,y(t),y(t-τ)) up to order Nβ exist and are continuous, then the local error of (56) admits the rigorous bound:(58)Li≤h(N+1)β1Γ((N+1)β)max⁡ϵ∈[0,1]y((N+1)β)(ts+ϵζih)∑j=1N‍cijmax⁡ε∈[0,1]fi((N)β)(εζih)+1Γ((N+1)β)∑j=1N‍cijmax⁡ϵ∈[0,1]fi((N)β)(ϵζih),and so (59)Li≤Ch(N+1)β.Proof.
Because of the order conditions and using (57), the first terms in the h expansion of (56) vanish. Then (58) is obtained and the proof is complete.Under the assumption of Theorem5, we see that the expressions (58) and (59) are bounded by a constant independent of h(N+1)β.Remark 6.
The proposed Runge-Kutta method with shifted Chebyshev polynomials (48) has order Nβ 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 thaty(t) is the exact solution of the system of FDDEs (20) and that ν(t) is the approximate solution. If the following inequalities are held:(60)ν(t0)-y(t0)≤ϕ(t0),t0∈[-τ,0],(61)Dt(β)ν(t)-f(t,ν(t),ν(t-τ))≤ε,(62)f(t,ν(t),ν(t-τ))-f(t,y(t),y(t-τ))≤lν(t)-y(t),then for t≥t0, one has the following error estimate: (63)νt-yt≤ϕ(t)exp⁡(l(tαΓ(α+1)-t0αΓ(α+1)))+εl(exp⁡(l(tαΓ(α+1)-t0αΓ(α+1)))-1).Proof.
For any chosen norm, we investigate the error estimated by(64)m(t)=ν(t)-y(t),and now, we try to estimate its natural growth.
Operate on both sides of (64) by the Caputo fractional derivative operator Dt(β) and using the triangle inequality yields(65)Dtβmt≤Dt(β)ν(t)-f(t,ν(t),ν(t-τ))+f(t,ν(t),ν(t-τ))-f(t,y(t),y(t-τ)),from Lipchitz condition (62), and we introduce δ(t)=Dt(β)ν(t)-f(t,ν(t),ν(t-τ)) which is called the defect of the approximate solution ν(t) such that δ(t)≤ε; we have (66)Dt(β)m(t)≤δ(t)+lm(t),δ(t)≤ε,(67)mt=ϕt,∀t∈[-τ,0].Now, we are about to solve instead of (67) the fractional differential equation and conclude that (68)m(t)≤u(t),(69)Dt(β)u(t)=lu(t)+ε,(70)ut=ϕt,∀t∈[-τ,0].To solve (69)-(70), we use the fractional complex transform T=tβ/Γ(β+1) and the modified Riemann-Liouville fractional derivative defined by(71)Dt(β)u(t)=1Γ(n-β)dndtn∫0t‍t-τn-β(u(τ)-u(0))dτ,where n-1≤β<n and n≥1.
Modified Riemann-Liouville fractional derivative has the following property:(72)Dt(β)u(T(t))=uT′(T)T(β)(t).Then (69)-(70) are transformed to(73)dudT=luT+ε,uT=ϕT,∀T∈[-τ,0].The exact solution for (73) can be written as (74)u(T)=ϕ(t)exp⁡(l(T-T0))+εl(exp⁡(l(T-T0))-1).Then the solution can be written using the original variables as follows:(75)ut=ϕ(t)exp⁡(l(tαΓ(α+1)-t0αΓ(α+1)))+εl(exp⁡(l(tαΓ(α+1)-t0αΓ(α+1)))-1).A combination of (68) and (75) yields the desired result (63).Theorem 8.
LetU be a neighborhood of (t,y(t)∣t0≤t≤tN) where y(t) is the exact solution of (20). Suppose that in U(76)dfdy≤l,μ(dfdy)≤land that the local error estimates (59) are valid in U. Then the global error E can be estimated by (77)E≤Γ(β)hNβtNβ-1·C′l(tNt01-βexp⁡(l(tNβΓ(β+1)-t0βΓ(β+1))-1)),where h=max⁡1<s≤Nhs-1,(78)C′=C,forl≥0;Cexp⁡(-lhβ),forl<0.And h is small enough for the numerical solution to remain in U.Proof.
Our task is to estimate the global error(79)E=y(tN)-yN,and this estimate is found along N-s steps of the numerical method (see Figure 1).
Estimate the local error using Theorem5 and Figure 1 to obtain (80)es≤Chs-1(N+1)β.From Theorem 7 with ε=0, we obtain (81)Es≤esexp⁡(l(tNβΓ(β+1)-tsβΓ(β+1))),from Figure 1, E=E1+E2+⋯+EN-1+EN,(82)E≤∑s=1N‍Es,as h=max⁡1<s≤Nhs-1, then we can deduce that (83)hs-1(N+1)β≤hNβhs-1β,and from (81)–(83), we deduce that(84)E≤ChNβ∑s=1N‍hs-1βexp⁡(l(tNβΓ(β+1)-tsβΓ(β+1)))=ChNβI,I≤∫t0tNexp⁡(l(tNβΓ(β+1)-tβΓ(β+1)))dt,forl≥0;∫t0tNexp⁡(l(tNβΓ(β+1)-hβ-tβΓ(β+1)))dt,forl<0.Then(85)I≤ΓβltNβ-1tNt01-β·exp⁡ltNβΓβ+1-t0βΓβ+1-1,forl≥0;Γ(β)ltNβ-1exp⁡(-lhβ)tNt01-β·exp⁡ltNβΓβ+1-t0βΓβ+1-1,forl<0.A combination of (84)-(85) yields the desired result (77).Figure 1
Global error estimation [42].Remark 9 (see [43]).
Since the global truncating errorE (77) tends to zero as h tends to zero, then the proposed method (43) is convergent; that is, |E|→0 as h→0.


## 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(86)Dtβyt=2Γ(3-β)t2-β-1Γ(2-β)t1-β+2τt-τ2-τ-y(t)+y(t-τ),y(t)=t2-t,t∈[-τ,0].The exact solution for this problem is y(t)=t2-t.For FDDEs (86) at (β=1, τ=1, h=0.25) and (β=0.1, τ=0.1, h=0.25) the absolute errors between the computed points and exact ones are presented in Figures 2 and 3 respectively.Figure 2
Plot of the absolute error withh=0.25, t∈[0,1], β=1, τ=1.Figure 3
Plot of the absolute error withh=0.25, t∈[0,1], β=0.1, τ=0.1exp⁡(-10t).Example 2.
Consider(87)Dt(3/4)y(t)=y(t-1)-y(t)+g(t),t>0,y(t)=sin(t),t∈[-1,0],where g(t) is defined by (88)gt=sint-sint-1+4t1/4Γ1/4F211;58,98;-t24,such that F21 is a special case of hypergeometric function.
The exact solution for this problem isy(t)=sin(t).
In  Table3, h=1/2n10, we present the results obtained with the proposed method at N=4 for different step sizes such that the maximum of the absolute error is given by (89)Error=max⁡j⁡yexacttj-ynumericaltj.
In the third column, the numerical convergence order which is defined bylog⁡2(Error(2h)/Error(h)) is introduced.Table 3
n
Error(h)
Numerical convergence order
0
0.001761
1
0.000264
2.73779
2
0.0000373
2.82329
3
4.821
×
1
0
-
6
2.95177
4
5.732
×
1
0
-
7
3.07222
5
6.850
×
1
0
-
8
3.06486Example 3.
Consider(90)Dt(1/2)y(t)=y(t-1)-y(t)+g(t),t>0,y(t)=sin(t),t∈[-1,0],where g(t) is defined by(91)gt=sint-sint-1+2tπF211;34,54;-t24.The exact solution for this problem is y(t)=sin(t).
In Table4, the computational results of the maximum of absolute errors are displayed at N=8 for different time steps h=1/2n10 with n=0,1,2,3,4,5. In the third column, the numerical convergence order is introduced also. These results are in concordance with our theoretical results.Table 4
n
Error(h)
Numerical convergence order
0
9.876
×
1
0
-
6
1
7.876
×
1
0
-
7
3.64839
2
5.866
×
1
0
-
8
3.74701
3
3.856
×
1
0
-
9
3.9272
4
1.846
×
1
0
-
10
4.38463
5
8.360
×
1
0
-
12
4.46476


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



---

*Source: 510875-2015-03-29.xml*