On the convergence of a new reliable algorithm for solving multi-order fractional differential equations

Accepted Manuscript

On the convergence of a new reliable algorithm for solving multi-order fractional differential equations Esmail Hesameddini, Azam Rahimi, Elham Asadollahifard PII: DOI: Reference:

S1007-5704(15)00351-2 10.1016/j.cnsns.2015.10.020 CNSNS 3678

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

17 August 2015 20 October 2015 23 October 2015

Please cite this article as: Esmail Hesameddini, Azam Rahimi, Elham Asadollahifard, On the convergence of a new reliable algorithm for solving multi-order fractional differential equations, Communications in Nonlinear Science and Numerical Simulation (2015), doi: 10.1016/j.cnsns.2015.10.020

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We introduce a new algorithm for solving multi-order fractional differential

CR IP T

equations. • A comparison between the presented method with some other well-known methods for solving M-FDEs is provided.

AC

CE

PT

ED

M

AN US

• We present an elegant way to show the convergence analysis of the RVIM.

1

ACCEPTED MANUSCRIPT

On the convergence of a new reliable algorithm for solving multi-order fractional differential equations✩

a Department

CR IP T

Esmail Hesameddinia,∗, Azam Rahimia , Elham Asadollahifarda of Mathematics, Shiraz University of Technology, P. O. Box 71555-313, Shiraz , Iran

Abstract

AN US

In this paper, we will introduce the reconstruction of variational iteration method (RVIM) to solve multi-order fractional differential equations (M-FDEs), which include linear and nonlinear ones. We will easily obtain approximate analytical solutions of M-FDEs by means of the RVIM based on the properties of fractional calculus. Moreover, the convergence of proposed method will be shown. Our scheme has been constructed for the fully general set of M-FDEs without

M

any special assumptions, and is easy to implement numerically. Therefore, our method is more practical and helpful for solving a broad class of M-FDEs. Nu-

ED

merical results are carried out to confirm the accuracy and efficiency of proposed method. Several numerical examples are presented in the format of table and graphs to make comparison with the results that previously obtained by some

PT

other well known methods.

Keywords: Multi-order fractional differential equation; Fractional calculus; Reconstruction of variational iteration method.

CE

2010 MSC: 34A08, 26A33, 65M12, 49M30

∗ Corresponding

AC

author Email addresses: [email protected] (Esmail Hesameddini), [email protected] (Azam Rahimi), [email protected] (Elham Asadollahifard)

Preprint submitted to Communications in nonlinear science and numerical simulationOctober 29, 2015

ACCEPTED MANUSCRIPT

1. Introduction The fractional calculus is a name for the theory of integrals and derivatives of arbitrary order. Fractional calculus generalizes the notion of derivative for those

5

it gives meaning to the expression

dα dtα f (t)

CR IP T

cases in which the differentiation order is not natural number; in other words,

in those cases where α is a fraction

or an irrational number. This generalization may be performed in several ways, leading to several sightly different definitions that do not always reach exactly

at the same results. The Riemann-Liouville operator and the Caputo operator commonly are used by authors. Over the last decades, the use of fractional

order derivatives has become more and more attractive in the broad field of

AN US

10

engineering to describe different kinds of problems. It is well known that the integer order differential operators are local, while the most important profit of using fractional differential equations (FDEs) is their nonlocal property. This means that the next state of a system depends not only upon its current state but also upon all of its historical states. Therefore, the memory effect of these

M

15

derivatives is one of the main reasons to use them in various applications. Since the fractional calculus is a powerful tool to describe physical systems that have

ED

a long-term memory, new possibilities appear in mathematics and theoretical physics. Therefore, FDEs have gained popularity for describing various phenomena, for instance in visco elasticity [1], colored noise [2], signal processing

PT

20

[2], control theory [3], anomalous diffusion [4] notably in chaotic systems [5] and in phase transitions [6]. In study of FDEs, one should note that finding an an-

CE

alytical or approximate solution is a challenging problem. Therefore, accurate methods for finding the solutions of FDEs are yet under investigation. Several

25

numerical methods for solving FDEs exist in the literature for example, Laplace

AC

transform method (Podlubny, 1999 [6]), Fourier transform method (Kemple and Beyer, 1997 [7]), Adomian’s decomposition method (Daftardar-Gejji and Jafari, 2005 [8] ; Daftardar-Gejji and Jafari, 2007 [9]), Homotopy analysis method (Liao, 2003 [10] ; Momani and Odibat, 2008 [11] ; Yildirim, 2009 [12] ; Ku-

30

mar and Singh, 2010 [13] ; Atangana and Secer, 2013 [14] ), Sumudu transform

3

ACCEPTED MANUSCRIPT

(Singh et al., 2011 [15] ; Sushila et al., 2013 [16] ; Atangana and Baleanu, 2013 [17] ), fractional Adams-Moulton method (Galeone and Garrappa, 2006 [18] ) and so on.

35

CR IP T

Multi-order fractional differential equations (M-FDEs) have been used to model various types of visco-elastic damping (see [6, 19, 20] ). Some numerical methods have been investigated for solving M-FDEs such as operational ma-

trix [21], Galerkin finite element method [22], predictor-corrector method [23], spectral collocation method [24] and Adams method [25]. Moreover, very few algorithms for the analytical solution of M-FDEs have been suggested [9, 26, 27] . And many of these methods are essentially used for particular types of these

AN US

40

equations, often just linear ones or even smaller classes. It should be noted that, most of these methods can not be generalized to nonlinear cases. The variational iteration method (VIM) that was first introduced by He [28] as a modification of the general Lagrange multiplier method [29] has been successfully applied to 45

many ordinary and partial differential equations [30, 31].

M

This encouraged researchers to extend this method for FDEs. So this extension is done and named as fractional variational iteration method (FVIM). Up

ED

to now, this method has been applied for solving different kinds of FDEs [32, 33]. Sweilam et al. have used FVIM for solving M-FDEs [34]. They converted the 50

equation into a system of FDEs and then applied FVIM to the resulting system.

PT

Not only these authors, but also some other authors did this strategy in order to solve M-FDEs. The shortcoming of this strategy is that, if the order of the

CE

M-FDE is very large, then the fractional differential system has many equations, so it is difficult to find analytical solution for this system. In 2015, we proposed a new alternative approach based on the variational

55

AC

iteration formulations and Laplace transform for solving fractional partial differential equations which was called the Reconstruction of Variational Iteration Method (RVIM) [35]. In the present work, we extend the RVIM for solving multi-order fractional

60

differential equations, which include linear and nonlinear ones. We can easily obtain approximate analytical solutions of M-FDEs by means of the RVIM and 4

ACCEPTED MANUSCRIPT

a few simple transformations which are based on the properties of the fractional calculus. Moreover, the convergence of proposed method will be studied. Our aim is to provide an accurate scheme that is robust, reliable, and reasonably inexpensive in terms of both set-up costs and the time taken to execute. Our

CR IP T

65

scheme has been constructed for the fully general set of M-FDEs without any special assumptions, and is easy to implement numerically. It is worth mention-

ing that the proposed method is capable of reducing the volume of the compu-

tational work as compared to some other classical methods. The outline of this 70

paper is as follows. Section 2 contains preliminaries. In section 3, the RVIM

AN US

is developed to solve the M-FDEs. Section 4 is devoted to the convergence of

proposed method. In section 5, extensive numerical experiments are presented to illustrate the accuracy and efficiency of our method. Finally, discussion and conclusion are summarized in section 6.

2. Preliminaries

M

75

In this section, we recall the basic definitions and operational properties of fractional integral and derivative. Many definitions of fractional calculus have

ED

been proposed in the past two centuries. These definitions include RiemannLiouville, Reize, Caputo and Grnwald-Letnikov fractional operators. The two most commonly used definitions are the Riemann-Liouville operator and the

PT

80

Caputo operator. In this part, we enlist some definitions and properties of the

CE

fractional calculus.

Definition 2.1. [36] A real function f (t), t > 0, is said to be in the space Cµ , µ ∈ R if there exists a real number p(> µ), such that f (t) = tp v(t), where v(t) ∈ C[0, ∞), and it is said to be in the space Cµm if f (m) ∈ Cµ ,m ∈ N .

AC

85

Definition 2.2. [36] The Riemann-Liouville fractional integral operator of order α ≥ 0, of a function f (t) ∈ Cµ , µ ≥ −1 is defined as: I α f (t) =

1 Γ(α)

Z

0

t

(t − τ )α−1 f (τ )dτ, α > 0, t > 0,

5

ACCEPTED MANUSCRIPT

such that I 0 f (t) = f (t). Note that Γ is the gamma function. For the Riemann-Liouville fractional integral we have the following properties Γ(β + 1) (t − a)β+α , Γ(β + α + 1)

I α I β f (t) = I β I α f (t) = I α+β f (t). 90

CR IP T

I α (t − a)β =

Definition 2.3. [36] The Riemann-Liouville fractional derivative of f , f (t) ∈ Cµ of order α ≥ 0 is defined as: α DR−L f (t) = Dn I n−α f (t).

AN US

Definition 2.4. [37] For a continuous function f (t), the Caputo fractional derivative of order α where n − 1 < α < n, is defined as Z t 1 (t − τ )n−α−1 Dn f (τ )dτ . Dcα f (t) = I n−α Dn f (t) = Γ(n − α) 0 It has the following property

M

I α Dcα f (t) = f (t) −

n X

k=0

f (k) (0+ )

xk . k!

Definition 2.5. The Laplace transform of the ordinary derivative Dn f (t) is given by

ED

95

PT

`{Dn f (t); s} = sn F (s) −

n−1 X

sn−k−1 f (k) (0),

k=0

where F (s) = `{f (t); s}.

CE

Furthermore, if f (k) (0) = 0, k = 0, 1, · · · , n − 1, then α 1)DR−L f (t) = Dcα f (t),

(1)

AC

2)Dcα I β f (t) = I β Dcα f (t) = Dcα−β f (t) = I β−α f (t),

(2)

and if n − 1 < α ≤ n, then Dcα tk = 0, k = 0, 1, · · · , n − 1.

The intrested readers should refer to [36, 37] for more properties of RiemannLiouville and Caputo fractional derivatives. 6

ACCEPTED MANUSCRIPT

100

3. Basic idea of the reconstruction of variational iteration method(RVIM) To implement the basic idea of our technique, we consider a general multi-

CR IP T

order fractional differential equation (M-FDE) as Dcαn x(t) = g(t, x(t), Dcα1 x(t), Dcα2 x(t), · · · , Dcαn−1 x(t)),

(3)

where 0 ≤ αi ≤ αn ≤ n, for i = 1, · · · , n − 1 and n − 1 < αn ≤ n, with the initial conditions

105

(4)

AN US

x(k) (0) = xk0 , k = 0, 1, · · · , n − 1,

and g : D = [0, T ] × R × R × · · · × R → R is a given continuous mapping and it has continuous and bounded partial derivatives γk = sup0≤t≤T |

∂g ∂vk

and

∂g |, k = 0, · · · , n − 1. ∂vk

M

Pn−1 xk If we take x(t) = x(t) − k=0 k!0 tk , then (3) can be written as  Pn−1 xk Pn−1 xk   Dcαn x(t) = g(t, x(t) + k=0 k!0 tk , Dcα1 x(t) + Dcα1 ( k=0 k!0 tk ), · · · ,    Pn−1 xk0 k  α α  Dc n−1 x(t) + Dc n−1 ( t )), k=0 k!

(6)

ED

       x(k) (0) = x(k) (0) − xk = 0, k = 0, 1, · · · , n − 1. 0

(5)

PT

According to (1) and the definition of Caputo derivative, we have Dcαn x(t) = αn DR−L x(t) = Dn I n−αn x(t). Let I n−αn x(t) = y(t), then problem (6) can be

written in the following form

CE

110

AC

y (n) (t) = g(t, Dcn−αn y(t) +

n−1 X k=0

n−1 X

+Dcα1 (

k=0

xk0 k k!

xk0 k n−αn +α1 t , Dc y(t) k! n−1 X

t ), · · · , Dcn−αn +αn−1 y(t) + Dcαn−1 (

k=0

xk0 k t )). k!

(7)

In variational iteration method(VIM), the correction functional is established by the general Lagrange multiplier which can be identified optimally via the variational theory. Not only RVIM requires no knowledge of variational theory but also without any restrictive assumptions we can derive an iterative relation. 7

ACCEPTED MANUSCRIPT

115

Indeed the RVIM is a new alternative approach based on the variational iteration formulations and the Laplace transform. To explain the RVIM, at first, by taking the Laplace transform from both sides of

sn `{y} = `{g(t, Dcn−αn y(t) +

n−1 X k=0

+Dcα1

n−1 X k=0

xk0 k

xk0 k n−αn +α1 t , Dc y(t) k!

t , · · · , Dcn−αn +αn−1 y(t) + Dcαn−1

k!

n−1 X k=0

xk0 k t )}. k!

(8)

AN US

So,

CR IP T

(7) and using the zero artificial initial conditions, the following result is obtained

n−1

`{y} =

X xk 1 0 k n−αn y(t) + `{g(t, D t , Dcn−αn +α1 y(t) c sn k! k=0

+Dcα1

n−1 X

xk0 k

k=0

Suppose that

1 sn

n−1 X k=0

xk0 k t )}. k!

(9)

= H(s). Then, by using the convolution theorem, we have

M

120

k!

t , · · · , Dcn−αn +αn−1 y(t) + Dcαn−1

`{y} = `{h ∗ g},

ED

where

`−1 {H(s)} = h(t) =

(10)

tn−1 . (n − 1)!

(11)

PT

Applying the inverse Laplace transform to both sides of (9), one obtains y(t) =

Z

h(t − τ )g(τ, Dcn−αn y(τ ) +

CE

0

t

+Dcα1

n−1 X k=0

xk0 k!

k

τ ,···

n−1 X k=0

xk0 k n−αn +α1 y(τ ) τ , Dc k!

, Dcn−αn +αn−1 y(τ )

+

Dcαn−1

n−1 X k=0

xk0 k τ )dτ. k!

(12)

AC

Thus, the following iteration formulation is resulted yi+1 (t) = y0 (t) +

Z

t

0

Dcn−αn +α1 yi (τ ) + Dcα1

h(t − n−1 X k=0

τ )g(τ, Dcn−αn yi (τ )

xk0 k!

+

n−1 X k=0

xk0 k τ , k!

τ k , · · · , Dcn−αn +αn−1 yi (τ ) + Dcαn−1

8

n−1 X k=0

xk0 k τ )dτ, k!

(13)

ACCEPTED MANUSCRIPT

where y0 (t) must be satisfied in the given initial conditions in order to impose 125

the actual conditions. Beginning with initial approximation y0 (t), the exact solution y(t) is obtained by y(t) = limi→∞ yi (t).

x(t) = Dcn−αn y(t) +

n−1 X k=0

xk0 k t . k!

CR IP T

Therefore, problem (3) has the following approximate solution

(14)

In the next section, we prove that the sequence {yi (t)}∞ i=1 , defined by (13) will converge to the solution of (3).

AN US

130

4. Convergence Analysis of RVIM for solving M-FDE

In this section, a convergence analysis will be given for the RVIM based on the error estimate. The main results are proposed in the following theorems. Theorem 4.1. Let y(t), yi (t) ∈ C n [0, T ], i = 0, 1, · · · . The error estimate is

M

135

ED

obtained by the following relation kEi+1 k∞ < kE0 k∞

(γnT αn )i+1 , Γ((i + 1)(αn − αn−1 ) + 1)

(15)

where γ = max0 ≤j≤n−1 γj , Ek (t) = yk (t) − y(t), k = 1, 2, · · · .

PT

Therefore, the sequence defined by (13) with y0 (t) = y0 converges to the exact solution of (3).

Proof. Evidently, from (13), we get

CE 140

AC

y(t) = y0 (t) + +Dcα1

n−1 X k=0

Z

t

0

xk0 k!

h(t − τ )g(τ, Dcn−αn y(τ ) +

n−1 X k=0

xk0 k n−αn +α1 y(τ ) τ , Dc k!

τ k , · · · , Dcn−αn +αn−1 y(τ ) + Dcαn−1

9

n−1 X k=0

xk0 k τ )dτ. k!

(16)

ACCEPTED MANUSCRIPT

Subtracting (13) from (16), results in Z

t

h(t − τ ){g(τ, Dcn−αn yi (τ ) +

0

+Dcα1

n−1 X k=0

xk0 k!

+Dcα1

k=0

xk0 k!

k=0

xk0 k n−αn +α1 τ , Dc yi (τ ) k!

τ k , · · · , Dcn−αn +αn−1 yi (τ ) + Dcαn−1

−g(τ, Dcn−αn y(τ ) + n−1 X

n−1 X

n−1 X

xk0 k!

k=0

τ k , Dcn−αn +α1 y(τ )

n−1 X k=0

n−1 X

τ k , · · · , Dcn−αn +αn−1 y(τ ) + Dcαn−1

xk0 k τ ) k!

CR IP T

Ei+1 (t) =

xk0 k τ )}dτ. k!

k=0

(17)

operator, one obtains Ei+1 (t) = I n {g(τ, Dcn−αn yi (τ ) +

k=0

xk0 k!

n−1 X k=0

xk0 k!

k=0

k!

n−1 X k=0

xk0 k τ ) k!

τ k , Dcn−αn +α1 y(τ )

ED

+Dcα1

xk0

k=0

xk0 k n−αn +α1 yi (τ ) τ , Dc k!

τ k , · · · , Dcn−αn +αn−1 yi (τ ) + Dcαn−1

−g(τ, Dcn−αn y(τ ) + n−1 X

n−1 X

M

+Dcα1

n−1 X

AN US

Using (11) and (17) and the definition of Riemann-Liouville fractional integral

τ k , · · · , Dcn−αn +αn−1 y(τ ) + Dcαn−1

n−1 X k=0

xk0 k τ )}. k!

(18)

tives

∂g ∂vk ,

k = 0, 1, · · · , n − 1. Using Lagrange’s theorem, we have

AC

CE

145

PT

Noting that g(t, v0 , v1 , · · · , vn−1 ) has continuous and bounded partial deriva-

10

ACCEPTED MANUSCRIPT

+Dcα1

n−1 X

xk0 k!

k=0

+Dcα1

xk0 k!

n−1 X k=0

xk0 k!

τ k , Dcn−αn +α1 y(τ )

τ k , · · · , Dcn−αn +αn−1 y(τ ) + Dcαn−1

k=0 [g20 (ξ(t))Dcn−αn Ei (t)

n−1 X

xk0 k τ ) k!

n−1 X

xk0 k τ )| k!

k=0

k=0

0 + · · · + gn+1 (ξ(t))Dcn−αn +αn−1 Ei (t)]|

AN US

= |I

n

k=0

xk0 k n−αn +α1 τ , Dc yi (τ ) k!

τ k , · · · , Dcn−αn +αn−1 yi (τ ) + Dcαn−1

−g(τ, Dcn−αn y(τ ) + n−1 X

n−1 X

CR IP T

|Ei+1 (t)| = |I n {g(τ, Dcn−αn yi (τ ) +

≤ γ0 |I n Dcn−αn Ei (t)| + γ1 |I n Dcn−αn +α1 Ei (t)| + · · · + γn−1 |I n Dcn−αn +αn−1 Ei (t)| ≤ γ0 I αn |Ei (t)| + γ1 I αn −α1 |Ei (t)| + · · · + γn−1 I αn −αn−1 |Ei (t)| = (γ0 I αn + γ1 I αn −α1 + · · · + γn−1 I αn −αn−1 )|Ei (t)| .. .

M

≤ (γ0 I αn + γ1 I αn −α1 + · · · + γn−1 I αn −αn−1 )i+1 |E0 (t)|

≤ (γ0 I αn + γ1 I αn −α1 + · · · + γn−1 I αn −αn−1 )i+1 max |E0 (s)| αn

αn −α1

αn −αn−1 i+1

0 ≤s≤T

+ ··· + I ) max |E0 (s)| 0 ≤s≤T Rt (t − τ )(i+1)αn −1 dτ ≤ γ i+1 max |E0 (s)|ni+1 0 Γ((i + 1)(αn − αn−1 )) 0 ≤s≤T (I

+I

ED

≤γ

i+1

PT

≤ γ i+1 max |E0 (s)|ni+1 0 ≤s≤T

t(i+1)αn , Γ((i + 1)(αn − αn−1 ))((i + 1)αn )

(19)

where gk0 is the partial derivative of function g for the k-th variable,

CE

and

AC

ξ(t) =

(t, Dcn−αn y(t)

+Dcα1

+

n−1 X k=0

n−1 X k=0

+Dcαn−1

xk0 k k!

n−1 X k=0

xk0 k t + θ(Dcn−αn Ei (t)), Dcn−αn +α1 y(t) k!

t + θ(Dcn−αn +α1 Ei (t)), · · · , Dcn−αn +αn−1 y(t)

xk0 k t + θ(Dcn−αn +αn−1 Ei (t))), 0 < θ < 1. k!

11

(20)

ACCEPTED MANUSCRIPT

Also

150

CR IP T

1 1 ≤ Γ((i + 1)(αn − αn−1 ))((i + 1)αn ) Γ((i + 1)(αn − αn−1 ))((i + 1)(αn − αn−1 )) 1 = , (21) Γ((i + 1)(αn − αn−1 ) + 1) and based on the convergence of Mittag-Leffler functions [6], we conclude that kEi+1 k∞ ≤ kE0 k∞

(γnT αn )i+1 . Γ((i + 1)(αn − αn−1 ) + 1)

(22)

Now since T, kE0 k∞ , γ and αn are constants and 0 < αn − αn−1 , if i → ∞ αn i+1

5. Numerical experiments

AN US

(γnT ) then kE0 k∞ Γ((i+1)(α → 0. This completes the proof. n −αn−1 )+1)

In this section, we apply our proposed method for solving some M-FDEs 155

to show its capability and efficiency. These examples are chosen because their

other well known methods.

M

closed form solutions are available or they have been solved previosly by some

Also, in order to show that our method is practical for two-dimensional M-FDEs,

160

ED

too, an example is presented.

Example 5.1. Consider the following nonlinear M-FDE

PT

D3 y(t) + Dc2.5 (y(t)) + y 2 (t) = t4 ,

(23)

CE

subject to the initial conditions y(0) = y 0 (0) = 0, y 00 (0) = 2.

(24)

The exact solution of this problem is y(t) = t2 [38, 39].

AC

Let y(t) = y(t) − t2 . Then, (23) can be rewrite as D3 y(t) + Dc2.5 (y(t)) + y 2 (t) + 2t2 y(t) = 0, y(0) = y 0 (0) = y 00 (0) = 0.

(25)

Applying the Laplace transform to (25), one obtains `{y(t)} =

1 `{−Dc2.5 (y(t)) − y 2 (t) − 2t2 y(t)}. s3 12

(26)

ACCEPTED MANUSCRIPT

Taking the inverse Laplace transform from both sides of (26), implies that Z t (t − τ )2 y(t) = (−Dc2.5 (y(τ )) − y 2 (τ ) − 2τ 2 y(τ ))dτ 2 0 = I 3 [−Dc2.5 (y(t)) − y 2 (t) − 2t2 y(t)]. (27) Therefore, the following iterative relation is obtained

CR IP T

165

y k+1 (t) = y 0 (t) + I 3 [−Dc2.5 (y k (t)) − y 2k (t) − 2t2 y k (t)],

(28)

where y 0 (t) = 0 and y k (t) indicates the k-th approximation of y(t).

y 0 (t) = 0,

AN US

According to (28), the following relations are obtained

y 1 (t) = I 3 [0] = 0, .. . y k (t) = 0.

170

M

Therefore, y(t) = limk→∞ y k (t) = 0. Thus, the closed form of solution of (23) is resulted as

ED

y(t) = y(t) + t2 = t2 , which is the exact solution of (23). In comparison with some other methods such as operational matrix of Chebyshev polynomials [38] and homotopy analysis

PT

method [39], we obtained this solution more easily.

CE

Example 5.2. Consider the following M-FDE 3.2t2.5 , Γ(0.5)

(29)

subject to the initial conditions

AC

175

D2 y(t) + Dc0.5 (y(t)) + y(t) = t3 + 6t +

y(0) = y 0 (0) = 0, t ∈ [0, 1].

(30)

The exact solution of this problem is y(t) = t3 [40]. Applying the Laplace transform to (29), one obtains `{y(t)} =

1 3.2t2.5 `{−Dc0.5 (y(t)) − y(t) + t3 + 6t + }. 2 s Γ(0.5) 13

(31)

ACCEPTED MANUSCRIPT

Table 1: Error estimation of RVIM for Example (5.2)

t

absolute error k=4

k=5

k=6

k=7

0

0

0

0

0.1

1.8e-11

0

0

0.2

3.771e-09

1.5e-11

0

0.3

8.625e-08

6.2e-10

0

0.4

8.0245e-07

9.03e-09

8.00e-11

0.5

4.5547e-06

7.27e-08

9.00e-10

0.6

188983e-05

4.019e-07

6.6e-09

1.00e-10

1.01e-12

0.7

6.31379e-05

1.7131e-06

3.60e-08

6.00e-10

7.03e-12

0.8

1.799522e-04

6.0333e-06

1.564e-07

3.2e-09

4.26e-11

0.9

4.541870e-04

1.83628e-05

5.735e-07

1.43e-08

2.18e-10

1

1.041373e-03

4.97998e-05

1.838e-06

5.44e-08

8.55e-10

CR IP T

k=3

0

0

0

0

0

0

0

0

0

0

0

M

AN US

0

PT

ED

Taking the inverse Laplace transform from both sides of (31), implies that Z t 3.2τ 2.5 y(t) = (t − τ )(−Dc0.5 (y(τ )) − y(τ ) + t3 + 6τ + )dτ Γ(0.5) 0 3.2t2.5 = I 2 [−Dc0.5 (y(t)) − y(t) + t3 + 6t + ]. (32) Γ(0.5) Therefore, the following iterative relation is obtained

CE

yk+1 (t) = y0 (t) + I 2 [−Dc0.5 (yk (t)) − yk (t) + t3 + 6t +

(33)

where y0 (t) = 0.

AC

180

3.2t2.5 ], Γ(0.5)

Table 1 shows the absolute errors of our method for different values of k and

t. These are obtained by using Maple software. We see that by increasing k, the numerical results become more and more accurate. In contrast with [40], where

185

the small step size h =

1 512

was used to get an absolute error of 8.14e − 07 at

the final time t = 1, we obtained the value of 8.55e − 10 for the absolute error. 14

ACCEPTED MANUSCRIPT

So, our method provide more accurate results. Example 5.3. Consider the following multi-order inhomogeneous Bagley-Torvik

CR IP T

fractional differential equation [41] D2 y(t) + Dc1.5 (y(t)) + y(t) = 1 + t, y(0) = 1, y 0 (0) = 1, t ∈ [0, 1]. 190

(34)

Let y(t) = y(t) − (1 + t). Then, (34) can be rewritten as follows

D2 y(t) + Dc1.5 (y(t)) + y(t) = 0, y(0) = 0, y 0 (0) = 0, t ∈ [0, 1].

`{y(t)} =

AN US

Applying the Laplace transform to (35), one obtains 1 `{−Dc1.5 (y(t)) − y(t)}. s2

(35)

(36)

Taking the inverse Laplace transform from both sides of (36), implies that Z t y(t) = (t − τ )(−Dc1.5 (y(τ )) − y(τ ))dτ = I 2 [−Dc1.5 (y(t)) − y(t)]. (37) 0

M

Therefore, the following iterative relation is obtained

y k+1 (t) = y 0 (t) + I 2 [−Dc1.5 (y k (t)) − y k (t)],

(38)

195

ED

where y 0 (t) = 0 and y k (t) indicates the k-th approximation for y(t). According to (38), the following relations are obtained

PT

y 0 (t) = 0,

CE

y 1 (t) = I 2 [0] = 0, .. . y k (t) = 0,

AC

therefore, y(t) = limk→∞ y k (t) = 0. Thus the closed form solution for (34) is resulted as y(t) = y(t) + 1 + t = 1 + t,

where this is the exact solution of (34). This example was also solved by the radial basis function (RBF) and approximate results were reported for various N 15

ACCEPTED MANUSCRIPT

200

[42]. We observe that for N = 12, in multiquadric radial basis function (MQ– RBF), Inverse MQ–RBF and Gaussian RBF, the absolute errors were obtained as 1.481e − 7, 4.682e − 5 and 5.123e − 8 respectively, whereas we obtained the

CR IP T

exact solution. Example 5.4. Finally, consider the following fractional partial differential 205

equation

∂u(x, t) ∂ 2γ u(x, t) , 0 ≤ γ ≤ 1, 0 ≤ x ≤ 1, =− ∂t ∂x2γ subject to the initial condition u(0, t) = et .

(39)

and cosγ (xγ ) = 210

Eγ (ixγ )+Eγ (−ixγ ) . 2

AN US

The exact solution of this problem is u(x, t) = et cosγ (xγ ) , where cosγ (xγ ) is the P∞ xkγ generalized cosine function defined by Mittag-function Eγ (xγ ) = k=0 Γ(kγ+1)

Let u(x, t) = u(x, t) − et . Then, (39) can be rewrite as follows

(40)

M

∂u(x, t) ∂ 2γ u(x, t) =− − et , u(0, t) = 0. ∂x2γ ∂t Applying the Laplace transform to (40), one obtains

ED

`{u(x, t)} =

1 ∂u(x, t) `{− − et }. s2γ ∂t

(41)

PT

Taking the inverse Laplace transform from both sides of (41), implies that Z x 1 ∂u(ζ, t) u(x, t) = (x − ζ)2γ−1 (− − et )dζ. (42) Γ(2γ) 0 ∂t

AC

CE

According to the RVIM, the recursive relation is given by  Rx 1   un+1 (x, t) = u0 (x, t) + Γ(2γ) (x − ζ)2γ−1 (− ∂un∂t(ζ,t) − et )dζ,  0      u (x, t) = 0. 0

16

(43)

ACCEPTED MANUSCRIPT

Then, we have x2γ , Γ(2γ + 1) x2γ x4γ u2 (x, t) = et (− + ), Γ(2γ + 1) Γ(4γ + 1) x2γ x4γ x6γ u3 (x, t) = et (− + − ), Γ(2γ + 1) Γ(4γ + 1) Γ(6γ + 1) .. . n X (−1)k x2kγ un (x, t) = et . Γ(2kγ + 1)

CR IP T

u1 (x, t) = −et

215

Therefore, u(x, t) =

∞ X

AN US

k=1

un (x, t) = et

k=1

∞ X (−1)k x2kγ . Γ(2kγ + 1)

(44)

k=1

Thus the closed form of solution for (39) is resulted as u(x, t) = u(x, t) + et = et

∞ X (−1)k x2kγ = et cosγ (xγ ). Γ(2kγ + 1)

M

k=0

Which is the exact solution of (39). This shows that this method can successfully be applied for solving multi-order fractional partial differential equations.

different values of γ are consistent with the exact solution of this equation.

PT

220

ED

Graphical results are depicted in Figures 1–3. These numerical solutions for

6. Conclusion

CE

In this paper, we have presented an alternative approach based on the VIM and laplace transform, named RVIM, to compute the solution of multi-order fractional differential equations (M-FDEs). The convergence of our method for solving M-FDEs has been discussed through a theorem. Not only RVIM requires

AC

225

no knowledge of variational theory in spite of VIM, but also without any restrictive assumptions could derive an iterative relation. From the computational point of view, the solutions obtained by our method were in excellent agreement with those obtained via previous works and also it was in very good conformity

230

with the exact solution. The obtained solution by using the suggested method 17

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

ED

M

Figure 1: Plot of the RVIM solution with γ = 0 and γ = 0.25 for Example (5.4).

AC

Figure 2: Plot of the RVIM solution with γ = 0.5 and γ = 0.75 for Example (5.4).

revealed that our method was more practical and helpful for solving a broad class of M-FDEs.

18

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 3: Plot of the RVIM solution with γ = 1 for Example (5.4).

7. Acknowledgments

The authors are indebted to the referees for their valuable comments and helpful suggestions which substantially improved the quality of this paper.

M

235

ED

References

[1] Flutter analysis of a flag of fractional viscoelastic material, Journal of Sound

PT

and Vibration 333.

[2] Fractional order control strategies for power electronic buck converters, Signal Process 86.

CE

240

[3] Y. W. T. Li, Y. Yang, Synchronization of fractional-order hyperchaotic systems via fractional-order controllers, Discrete Dyn. Nat. Soc. 2014 (2014)

AC

1–14. doi:10.1016/0003-4916(63)90068-X.

[4] P. Cosenza, D. Koroak, Secondary consolidation of clay as an anomalous

245

diffusion process, International Journal for Numerical and Analytical Methods in Geomechanics 38 (12) (2014) 1231–46. doi:10.1002/nag.2256/ abstract.

19

ACCEPTED MANUSCRIPT

[5] R. Zhang, J. Gong, Synchronization of the fractional-order chaotic system via adaptive observer, Systems Science and Control Engineering 2 (2014) 751–54. doi:10.1080/21642583.2014.891955.VdHU5rkzSP8.

250

CR IP T

[6] Fractional differential equations, Academic Press, New York.

[7] S. Kemple, H. Beyer, Global and causal solutions of fractional differential

equations. in: Transform methods and special functions: Varna96, Proceedings of 2nd International Workshop (SCTP), Singapore 19 (1–12) (1997) 210–16. doi:10.1080/00207160701405477.

255

AN US

[8] V. Daftardar-Gejji, H. Jafari, Adomian decomposition: A tool for solving a

system of fractional differential equations, J. Math. Anal. Appl. 301 (2005) 508–18. doi:10.1016/j.jmaa.2004.07.039.

[9] Solving a multi-order fractional differential equation using adomian decomposition, J. Math. Anal. Appl. 189.

260

M

[10] Beyond perturbation: Introduction to the homotopy analysis method, CRC

ED

Press/Chapman and Hall, Boca Raton. [11] S. Momani, Z. Odibat, A novel method for nonlinear fractional partial differential equations: Combination of dtm and generalized taylor’s formula, J. Comput. Appl. Math. 220 (1–2) (2008) 85–95. doi:10.1016/j.cam.

PT

265

2007.07.033.

CE

[12] A. Yildirim, An algorithm for solving the fractional nonlinear schrondinger equation by means of the homotopy perturbation method, Int. J. Nonlinear Sci. Num. Simulat. 10 (2009) 445–51. doi:10.1515/IJNSNS.2009.10.4.

AC

270

445.

[13] A fractional model of gas dynamics equation by using laplace transform, Zeitschrift fur Naturforschung 67a.

[14] The time-fractional coupled-kortewegde-vries equations, Abst. Appl. Anal.

20

ACCEPTED MANUSCRIPT

[15] Homotopy perturbation sumudu transform method for nonlinear equations, Adv. Theor. Appl. Mech. 4.

275

[16] J. S. Sushila, Y. Shihodia, An efficient analytical approach for mhd viscous

CR IP T

flow over a stretching sheet via homotopy perturbation sumudu transform method, Ain Shams Eng. J. 4 (2013) 549–55. doi:10.1016/j.asej.2012. 12.002. 280

[17] Nonlinear fractional jaulent-miodek and whitham-broer-kaup equations within sumudu transform, Abstract Appl. Anal.

AN US

[18] On multistep methods for differential equations of fractional order, Mediterr. J. Math. 3.

[19] R. Bagley, P. Torvik, On the appearance of the fractional derivative in the behaviour of real materials, J. Appl. Mech. 51 (1984) 294–8. doi:

285

10.1115/1.3167615.

M

[20] Finite element modelling of viscoelastic materials on the theory of fractional

ED

calculus, Ph.D. Thesis, Pennsylvania State University. [21] M. A. D. Rostamy, H. Jafari, C. Khalique, Computational method based on bernstein operational matrices for multi-order fractional differential equa-

290

PT

tions, Filomat 3 (1) (2014) 591–601. doi:10.2298/FIL1403591R. [22] The galerkin finite element method for a multi-term time-fractional diffu-

CE

sion equation, J. Comput. Phys. 281. [23] N. F. K. Diethelm, A. D. Freed, A predictor-corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynam.

AC

295

29 (2002) 3–22. doi:10.1023/A:1016592219341.

[24] F. Ghoreishi, P. Mokhtary, Spectral collocation method for multi-order fractional differential equations, International Journal of Computational Methods 55 (10) (2014) 23 pages. doi:10.1142/S0219876213500722.

21

ACCEPTED MANUSCRIPT

300

[25] K. Diethelm, N. Ford, Multi-order fractional differential equations and their numerical solution, Appl. Math. Comput. 154 (3) (2004) 621–40. doi: 10.1016/S0096-3003(03)00739-2.

CR IP T

[26] I. T. H. Jiang, F. Liu, K. Burrage, Analytical solutions for the multi-term time-space caputo-riesz fractional advection-diffusion equations on a finite

domain, J. Math. Anal. Appl. 389 (2) (2012) 1117–27. doi:10.1016/j.

305

jmaa.2011.12.055.

[27] Analytical solutions for the multi-term time-space fractional advection-

AN US

diffusion equations with mixed boundary conditions, Nonlinear Anal. Real World Appl. 14. 310

[28] Some applications of nonlinear fractional differential equations and their approximations, Bull. Sci. Technol. 15 (12) (1999) 86–90.

[29] General use of the lagrange multiplier in non-linear mathematical physics,

M

in: Variational methods in the mechanics of solids, Pergamon Press, New York.

[30] F. K. M.T. Darvishi, A. Soliman, The numerical simulation for stiff systems

ED

315

of ordinary differential equations, Comput. Math. Appl. 54 (2007) 1055–63. doi:10.1016/j.camwa.2006.12.072.

PT

[31] G. Draganescu, Application of a variational iteration method to linear and nonlinear viscoelastic models with fractional derivatives, J. Math. Phys. 47 (8) (2006) 802–902. doi:10.1063/1.2234273.

CE

320

[32] G. Wu, D. Baleanu, Variational iteration method for the burgers flow with

AC

fractional derivatives-new lagrange multipliers, Appl. Math. Model. 37 (9) (2013) 6183–90. doi:10.1016/j.apm.2012.12.018.

[33] H. T. H. Jafari, D. Baleanu, A modified variational iteration method

325

for solving fractional riccati differential equation by adomian polynomials, Fract. Calc. Appl. Anal. 16 (1) (2013) 109–22. s13540-013-0008-9. 22

doi:10.2478/

ACCEPTED MANUSCRIPT

[34] M. K. N.H. Sweilam, R. Al-Bar, Numerical studies for a multi-order fractional differential equation, Phys. Lett. A 371 (2007) 26–33. doi: 10.1016/j.physleta.2007.06.016.

330

CR IP T

[35] E. Hesameddini, A. Rahimi, Solving fractional partial differential equations with variable coefficients by the reconstruction of variational iteration method, Z. Naturforsch. 5a (70) (2015) 375–82. doi:10.1515/ zna-2015-1017. 335

[36] Theory and applications of fractional differential equations, San Diego: El-

AN US

sevier.

[37] The fractional calculus: theory and application of differentiation and integration to arbitrary order, Academic Press.

[38] M. Seifollahi, A. Shamloo, Numerical solution of nonlinear multi-order fractional differential equations by operational matrix of chebyshev polynomi-

340

M

als, World Appl. Program. 3 (3) (2013) 85–92. doi:10.1155/2012/510180. [39] S. D. H. Jafari, H. Tajadodi, Solving a multi-order fractional differential

ED

equation using homotopy analysis method, J. King Saud Univ. Sci. 23 (2011) 151–55. doi:10.1016/j.jksus.2010.06.023. [40] N. Ford, A. Joseph, Connolly,systems-based decomposition schemes for

PT

345

the approximate solution of multi-term fractional differential equations, J. Comput. Appl. Math. 229 (2009) 382–91. doi:10.1016/j.cam.2008.

CE

04.003.

AC

[41] Numerical solution of the bagley-torvik equation, BIT 42.

350

[42] F. G. B. Fakhr Kazemi, Error estimate in fractional differential equations using multiquadratic radial basis functions, J. Comput. Appl. Math. 245 (2013) 133–47. doi:10.1016/j.cam.2012.12.011.

23