Haar wavelet method for solving fractional partial differential equations numerically

Applied Mathematics and Computation 227 (2014) 66–76

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Haar wavelet method for solving fractional partial differential equations numerically Lifeng Wang, Yunpeng Ma ⇑, Zhijun Meng School of Aeronautic Science and Technology, Beihang University, Beijing 100191, China

a r t i c l e

i n f o

a b s t r a c t

Keywords: Haar wavelet Operational matrix Fractional partial differential equations Sylvester equation Numerical solution

In this paper, a wavelet operational method based on Haar wavelet is proposed to solve the fractional partial differential equations in the Caputo derivative sense. We give the Haar wavelet operational matrix of fractional order integration. A truncated Haar wavelet series together with the wavelet operational matrix are utilized to reduce the fractional partial differential equations to Sylvester equations. In addition, some examples are presented to show the efficiency and the accuracy of the approach. Crown Copyright Ó 2013 Published by Elsevier Inc. All rights reserved.

1. Introduction Fractional calculus is an old mathematical topic with history as long as integer order calculus. However, its practical applications have recently been investigated [1–3]. Many physics and engineering systems can be elegantly modeled with the help of the fractional derivative, such as dielectric polarization [4], electrolyte–electrolyte polarization [5], and viscoelastic systems [6]. Owing to the increasing applications, there has been important interest in developing numerical methods for the solution of fractional differential equations. These methods include adomian decomposition method [7,8], generalized differential transform method [9,10], variational iteration method [11], finite difference method [12–13], homotopy analysis method [14], and wavelet method [15–17]. A small number of algorithms for the numerical solution of fractional partial differential equations have been proposed in recent years. Podlubny [18] used the Laplace transform method to solve the fractional partial differential equations with constant coefficients. Zaid Odibat and Shaher Momani [19] applied generalized differential transform method to solve the numerical solution of linear partial differential equations of fractional order. Zhang [20] discussed a practical implicit method to solve a class of initial boundary value space–time fractional convection–diffusion equations with variable coefficients. In this paper, we consider a class of fractional partial differential equations

@au @bu þ ¼ f ðx; tÞ; @xa @tb

ð1Þ

subject to

 @u @t 

¼ u1 ðtÞ;

x¼0

uð0; tÞ ¼ h1 ðtÞ;

 @u @x 

¼ u2 ðxÞ;

ð2Þ

t¼0

uðx; 0Þ ¼ h2 ðxÞ;

⇑ Corresponding author. E-mail addresses: [email protected] (L. Wang), [email protected], [email protected] (Y. Ma), [email protected] (Z. Meng). 0096-3003/$ - see front matter Crown Copyright Ó 2013 Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.11.004

ð3Þ

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76 a

67

b

uðx;tÞ uðx;tÞ where @ @x and @ @t are fractional derivative of Caputo sense, f ; u1 ; u2 ; h1 ; h2 are the known continuous functions, uðx; tÞ is a b the unknown function, 0 < a; b 6 1. The paper is organized as follows: In Section 2 we introduce some necessary definitions and mathematical preliminaries of fractional calculus. In Section 3, after describing the basic formulation of Haar wavelet, we give the Haar wavelet operational matrix of fractional integration. In Section 4 we provide several examples to show the efficiency and simplicity of the method. Concluding remarks are given in the last section.

2. Preliminaries and notations In this section, we present some definitions, notations and preliminaries of the fractional calculus theory which will be used in this work [18]. Definition 2.1. The Riemann–Liouville fractional integral operator of order a P 0 of a function is defined as

J a uðtÞ ¼

1 CðaÞ

Z

t

ðt  sÞa1 uðsÞds;

a > 0; t > 0;

ð4Þ

0

J 0 uðtÞ ¼ uðtÞ:

ð5Þ a

The properties of the operator J are given as follows

(i) J a J b uðtÞ ¼ J aþb uðtÞ, (ii) J a J b uðtÞ ¼ J b J a uðtÞ, aþc (iii) J a tc ¼ CCðaðþcþ1Þ . cþ1Þ t

Definition 2.2. The fractional derivative of uðtÞ in the Caputo sense is defined as a

D uðtÞ ¼

8 r < d uðtÞ r ; dt

:

1

CðraÞ

a ¼ r 2 N; Rt

uðrÞ ðsÞ 0 ðtsÞarþ1

ds; 0 6 r  1 < a < r:

ð6Þ

The Caputo fractional derivative of order a is also defined as Da uðtÞ ¼ J ra Dr uðtÞ, where Dr is the usual integer differential operator of order r.

3. Haar wavelet and operational matrix of the fractional integration 3.1. Haar wavelet For t 2 ½0; 1, the orthogonal set of Haar wavelet functions are defined by [21]:

1 h0 ðtÞ ¼ pffiffiffiffiffi ; m

ð7Þ

8 j=2 k1 > 6 t < k1=2 >2 ; 2j 2j 1 < k hi ðtÞ ¼ pffiffiffiffiffi 2j=2 ; k1=2 6 t < 2j 2j m> > : 0; otherwise;

ð8Þ

where i ¼ 0; 1; 2; . . . ; m  1, m ¼ 2pþ1 and p is a positive integer. j and k represent the integer decomposition of the index i, i.e. i ¼ 2j þ k  1. 3.2. Function approximation A function uðtÞ 2 L2 ð½0; 1ÞÞ can be expanded into Haar wavelet by

uðtÞ ¼

1 X ci hi ðtÞ; i¼0

R1 where ci ¼ 0 uðtÞhi ðtÞdt. If uðtÞ is approximated as piecewise constant during each subinterval, Eq. (9) will be terminated at finite terms

ð9Þ

68

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

uðtÞ 

m1 X

ci hi ðtÞ ¼ cT Hm ðtÞ;

ð10Þ

i¼0 T

where c ¼ ½c0 ; c1 ; . . . ; cm1 T ; Hm ðtÞ ¼ ½h0 ðtÞ; h1 ðtÞ; . . . ; hm1 ðtÞ , m is a power of 2. The matrix form of Eq. (10) is

u ¼ cT H;

ð11Þ pþ1

where the row vector u is the discrete form of the function uðtÞ. H is Haar wavelet matrix of order m ¼ 2 i.e.

h0 ðt 0 Þ

h0 ðt 1 Þ



h0 ðt m1 Þ

3

6 h ðt Þ 6 1 0 H¼6 .. 6 4 .

h1 ðt1 Þ .. .

 .. .

h1 ðt m1 Þ .. .

7 7 7: 7 5

2

; p ¼ 0; 1; 2; . . .,

hm1 ðt0 Þ hm1 ðt1 Þ    hm1 ðtm1 Þ For arbitrary function uðx; tÞ 2 L2 ð½0; 1Þ  ½0; 1ÞÞ, it can be also expanded into Haar series by [15]

uðx; tÞ 

m1X m1 X

uij hi ðxÞhj ðtÞ;

ð12Þ

i¼0 j¼0

where uij ¼ hhi ðxÞ; huðx; tÞ; hj ðtÞii, hhi ðxÞ; hj ðxÞi ¼ Eq. (12) will be written as

R1 0

hi ðxÞhj ðxÞdx.

uðx; tÞ  HTm ðxÞUHm ðtÞ:

ð13Þ

In this paper, we apply wavelet collocation method to determine the coefficients uij . These collocation points are shown in the following

xl ¼ t l ¼ ðl  1=2Þ=m;

l ¼ 1; 2; . . . ; m:

ð14Þ

Discreting Eq. (13) by the step Eq. (14), we obtain the matrix form of Eq. (13)

C ¼ HT UH;

ð15Þ

where U ¼ ½uij mm , C ¼ ½uðxi ; tj Þmm . From the definition of Haar wavelet functions, we may know that H is a orthogonal matrix. 3.3. Convergence of the Haar wavelet bases In this part, we assume that

9M > 0;

@uðx;tÞ @x

is continuous and bounded on ð0; 1Þ  ð0; 1Þ, then

  @uðx; tÞ    @x  6 M:

8x; t 2 ð0; 1Þ  ð0; 1Þ;

ð16Þ

Suppose um ðx; tÞ is the following approximation of uðx; tÞ,

um ðx; tÞ ¼

m1 Xm1 X

unl hn ðxÞhl ðtÞ:

ð17Þ

n¼0 l¼0

Then we have

uðx; tÞ  um ðx; tÞ ¼

1 X 1 1 1 X X X unl hn ðxÞhl ðtÞ ¼ unl hn ðxÞhl ðtÞ: n¼m l¼m

ð18Þ

n¼2pþ1 l¼2pþ1

Theorem 1. Suppose that the functions um ðx; tÞ obtained by using Haar wavelet are the approximation of uðx; tÞ, then we have the error bounded as follows

M 1 kuðx; tÞ  um ðx; tÞkE 6 pffiffiffi 3 ; 3m where kuðx; tÞkE ¼ ð

R1 R1 0

0

u2 ðx; tÞdxdtÞ

1=2

.

Proof. The orthonormality of the sequence fhi ðxÞg on ½0; 1Þ implies that

69

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

Z



1

hn ðxÞhn0 ðxÞdx ¼

1=m; n ¼ n0

ð19Þ

n–n0 :

0;

0

Therefore

kuðx; tÞ  um ðx; tÞk2E ¼ ¼

Z

1

Z

1

½uðx; tÞ  um ðx; tÞ2 dxdt

0 0 1 1 X X

1 1 X X

unl un0 l0

Z

n¼2pþ1 l¼2pþ1 n0 ¼2pþ1 l0 ¼2pþ1

1

Z hn ðxÞhn0 ðxÞdx

0

0

1

hl ðtÞhl0 ðtÞdt

 ¼

1 1 X 1 X u2nl ; 2 m pþ1 pþ1 n¼2

ð20Þ

l¼2

where unl ¼ hhn ðxÞ; huðx; tÞ; hl ðtÞii. According to Eqs. (7) and (8), we have

huðx; tÞ; hl ðtÞi ¼

Z 0

1

2j=2 uðx; tÞhl ðtÞdt ¼ pffiffiffiffiffi m

Z

ðk12Þ2j

ðk1Þ2j

uðx; tÞdt 

Z

!

k2j

uðx; tÞdt : ðk12Þ2j

ð21Þ

Using mean value theorem of integrals:

    1 1  2j ; k   2j 6 t 2 < k  2j ; 9t 1 ; t 2 : ðk  1Þ  2j 6 t1 < k  2 2 such that

2j=2 huðx; tÞ; hl ðtÞi ¼ pffiffiffiffiffi m

 k

      1 j 1 j 2  ðk  1Þ2j uðx; t 1 Þ  k2j  k  2 uðx; t2 Þ 2 2

2j=21 ¼ pffiffiffiffiffi ðuðx; t 1 Þ  uðx; t2 ÞÞ; m

ð22Þ

hence

*

+ Z 2j=21 2j=21 1 hn ðxÞ; pffiffiffiffiffi ðuðx; t 1 Þ  uðx; t2 ÞÞ ¼ pffiffiffiffiffi hn ðxÞðuðx; t1 Þ  uðx; t 2 ÞÞdx m m 0 Z 1  Z 1 2j=21 hn ðxÞuðx; t 1 Þdx  hn ðxÞuðx; t2 Þdx ¼ pffiffiffiffiffi m 0 0 ! Z ðk1Þ2j Z k2j Z ðk1Þ2j Z k2j 2 2 1 ¼ uðx; t 1 Þdx  uðx; t1 Þdx  uðx; t 2 Þdx þ uðx; t 2 Þdx : 2m ðk1Þ2j ðk12Þ2j ðk1Þ2j ðk12Þ2j

unl ¼

Using mean value theorem of integrals again:

    1 1  2j ; k   2j 6 x2 ; x4 < k  2j ; 9x1 ; x2 ; x3 ; x4 : ðk  1Þ  2j 6 x1 ; x3 < k  2 2 such that

          1 1 j 1 j 1 j 2  ðk  1Þ2j uðx1 ; t 1 Þ  k2j  k  2 uðx2 ; t1 Þ  k  2  ðk  1Þ2j uðx3 ; t 2 Þ k 2m 2 2 2     1 j 1 j 2 uðx4 ; t 2 Þ ¼ jþ2 ½ðuðx1 ; t1 Þ  uðx2 ; t1 ÞÞ  ðuðx3 ; t 2 Þ  uðx4 ; t 2 ÞÞ; ð23Þ þ k2  k  2 2 m

unl ¼

hence

u2nl ¼

1 2

2jþ4

m2

½ðuðx1 ; t 1 Þ  uðx2 ; t 1 ÞÞ  ðuðx3 ; t2 Þ  uðx4 ; t2 ÞÞ2 :

ð24Þ

Using mean value theorem of derivatives:

9n1 ; n2 : x1 6 n1 < x2 ; x3 6 n2 < x4 ; such that

2 @uðn1 ; t 1 Þ @uðn2 ; t2 Þ ðx  x Þ  ðx  x Þ 2 1 4 3 @x @x 22jþ4 m2 (   )  2  2 @uðn1 ; t 1 Þ@uðn2 ; t 2 Þ 1 2 @uðn1 ; t 1 Þ 2 @uðn2 ; t 2 Þ   : 6 2jþ4 ðx2  x1 Þ þ ðx4  x3 Þ þ 2ðx2  x1 Þðx4  x3 Þ   @x @x @x @x 2 m2

u2nl ¼

1



ð25Þ

70

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

Combining Eqs. (16) and (25), we obtain

4M2

u2nl 6

2

4jþ4

m2

¼

M2 2

4jþ2

m2

ð26Þ

:

Substituting Eq. (26) into Eq. (20), then we have

kuðx; tÞ 

um ðx; tÞk2E

0 1 0 1 jþ1 jþ1 jþ1 jþ1 2 1 1 1 2X 12X 1 1 2X 12X 1 X X 1 X 1 X 1 M 2 2 @ @ A ¼ 2 unl ¼ 2 unl A 6 2 m m j¼pþ1 m j¼pþ1 24jþ2 m2 n¼2pþ1 l¼2pþ1 n¼2j l¼2j n¼2j l¼2j 0 1 jþ1 jþ1 2 1 2X 12X 1 M2 X 1 1 M2 1 @ A¼ M ¼ : ¼ 4 2ðpþ1Þ 4jþ2 4 m j¼pþ1 3m 2 3 m6 2 n¼2j l¼2j

ð27Þ

Therefore

M 1 kuðx; tÞ  um ðx; tÞkE 6 pffiffiffi 3 : 3m

ð28Þ

From the Eq. (28), we can find that kuðx; tÞ  um ðx; tÞkE ! 0 when m ! 1. The larger the value of m, the more accurate the numerical solution. h 3.4. Operational matrix of the fractional integration In this part, we may simply introduce the operational matrix of fractional integration of Haar wavelet, more detailed introduction can be found in Ref. [22]. If J a is fractional integration operator of Haar wavelet, we get:

ðJ a Hm ÞðtÞ ¼ Pamm Hm ;

ð29Þ

a

where Pmm is called operational matrix of fractional integration of Haar wavelet. We can also know from Ref. [22]:

Pamm ¼ HF amm HT ;

ð30Þ

where

2 a

F mm

1

n1

n2

6 1 n1 6 0 6 1 1 6 0 ¼ a 0 1 m Cða þ 2Þ 6 6         4 0 0 0 aþ1

aþ1

   nm1

3

7    nm2 7 7    nm3 7 7; 7   5  1

ð31Þ

aþ1

and nk ¼ ðk þ 1Þ  2k þ ðk  1Þ . For instance, if a ¼ 0:5; m ¼ 8, we have

3 0:7523 0:2203 0:1558 0:0820 0:1102 0:0580 0:0447 0:0377 6 0:2203 0:3116 0:1558 0:2296 0:1102 0:0580 0:1756 0:0782 7 7 6 7 6 6 0:0410 0:1148 0:2203 0:0350 0:1102 0:1623 0:0389 0:0063 7 7 6 6 0:0779 0:0779 0 0:2203 0 0 0:1102 0:1623 7 7 6 ¼6 7: 7 6 0:0094 0:0196 0:0812 0:0032 0:1558 0:0247 0:0026 0:0009 7 6 6 0:0112 0:0439 0:0551 0:0194 0 0:1558 0:0247 0:0026 7 7 6 7 6 4 0:0145 0:0145 0 0:0812 0 0 0:1558 0:0247 5 2

P0:5 88

0:0275 0:0275

0

0:0551

0

0

0

0:1558

4. Applications and results In this section, we will use the Haar wavelet operational matrix of fractional order integration to solve the fractional partial differential equation Eq. (1). To demonstrate the effectiveness of this method, we consider four numerical examples. Example 1. Consider the following nonhomogeneous partial differential equation

@ 1=4 u @ 1=4 u þ ¼ f ðx; tÞ; @x1=4 @t1=4

0 6 x; t 6 1;

ð32Þ

71

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

  3=4 tþxt 3=4 Þ   such that @u ¼ @u ¼ uð0; tÞ ¼ uðx; 0Þ ¼ 0, f ðx; tÞ ¼ 4ðx3Cð3=4Þ , the exact solution is xt. @x t¼0 @t x¼0 T @2 u Let @x@t  Hm ðxÞUHm ðtÞ, then

@u ¼ @x @u ¼ @t

Z

t

0

Z

x

0

  Z t @2u @u @u dt þ   ½HTm ðxÞUHm ðtÞdt þ  ¼ HTm ðxÞUP1mm Hm ðtÞ; @x@t @x t¼0 @x t¼0 0

ð33Þ

  Z x T @2u @u @u  ½HTm ðxÞUHm ðtÞdx þ  ¼ HTm ðxÞ½P1mm  UHm ðtÞ: dx þ  @x@t @t x¼0 @t x¼0 0

ð34Þ

Therefore T

T

uðx; tÞ  HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ uð0; tÞ ¼ HTm ðxÞ½P 1mm  UP 1mm Hm ðtÞ:

ð35Þ

Then we have

  T @ 1=4 u 3=4 @u 1  J 3=4 ðHTm ðxÞUP1mm Hm ðtÞÞ ¼ HTm ðxÞ½P3=4 ¼ J mm  UP mm H m ðtÞ; @x1=4 @x

ð36Þ

  T T @ 1=4 u 3=4 @u 3=4  J 3=4 ðHTm ðxÞ½P1mm  UHm ðtÞÞ ¼ HTm ðxÞ½P1mm  UPmm ¼ J Hm ðtÞ: @t @t1=4

ð37Þ

Similarly, f ðx; tÞ may be expanded by the Haar wavelet functions as follows

f ðx; tÞ  HTm ðxÞFHm ðtÞ;

ð38Þ

where F ¼ ½fij mm .Substituting Eqs. (36) and (37) and Eq. (38) into Eq. (32), we have T

T

1 T 1 3=4 T HTm ðxÞ½P3=4 mm  UP mm Hm ðtÞ þ H m ðxÞ½P mm  UP mm H m ðtÞ ¼ H m ðxÞFH m ðtÞ:

ð39Þ

Dispersing Eq. (39) by the points ðxi ; t j Þ, i ¼ 1; 2; . . . ; m and j ¼ 1; 2; . . . ; m, we obtain T

T

1 T 1 3=4 T HT ½P3=4 mm  UP mm H þ H ½P mm  UP mm H ¼ H FH:

ð40Þ

Eq. (40) can be also written as T

T

T

3=4 3=4 1 1 1 ½P1 mm  ½P mm  U þ UP mm P mm ¼ ½P mm  FP mm :

ð41Þ

Eq. (41) is a Sylvester equation which is solved easily by using Matlab software. Solving it, we can get U. Then using Eq. (35), we obtain the approximation uðx; tÞ. The numerical results for m ¼ 8, m ¼ 16, m ¼ 32 are shown in Figs. 1–3. The exact solution is shown in Fig. 4. From the Figs. 1–4, we may see clearly that numerical solutions are in very good agreement with exact solution. Example 2. Consider the following fractional partial differential equation

0 6 x; t 6 1;

ð42Þ

1 0.8 0.6 u(x,t)

@ 1=3 u @ 1=2 u þ ¼ f ðx; tÞ; @x1=3 @t1=2

0.4 0.2 0 1 0.5

t

0

0

0.2

0.4

x

Fig. 1. Numerical solution of m ¼ 8.

0.6

0.8

1

72

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

subject to



@u @x t¼0

¼ 2x;



@u @t x¼0

5=3

3=2

uðx; 0Þ ¼ x2 , f ðx; tÞ ¼ CCð3Þx þ CCð3Þt . ð8=3Þ ð5=2Þ

uð0; tÞ ¼ t 2 ;

¼ 2t;

The exact solution of the equation is x2 þ t 2 . Let

@2 u @x@t

 HTm ðxÞUHm ðtÞ, then

@u ¼ @x @u ¼ @t Hence

Z

t

0

Z

x

0

  Z t @2u @u @u T dt þ   ½Hm ðxÞUHm ðtÞdt þ  ¼ HTm ðxÞUP1mm Hm ðtÞ þ 2x; @x@t @x t¼0 @x t¼0 0

ð43Þ

  Z xh i h iT @2u @u @u dx þ   HTm ðxÞUHm ðtÞ dx þ  ¼ HTm ðxÞ P1mm UHm ðtÞ þ 2t: @x@t @t x¼0 @t 0 x¼0

ð44Þ

h iT T uðx; tÞ  HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ x2 þ uð0; tÞ ¼ HTm ðxÞ P1mm UP1mm Hm ðtÞ þ x2 þ t 2 :

ð45Þ

Then we have

  h iT @ 1=3 u 2Cð2Þ 5=3 2=3 @u 1  J 2=3 ðHTm ðxÞUP1mm Hm ðtÞ þ 2xÞ ¼ HTm ðxÞ P2=3 ¼ J x ; mm UP mm H m ðtÞ þ 1=3 @x @x Cð8=3Þ

ð46Þ

  T T @ 1=2 u @u 2Cð2Þ 3=2  J 1=2 ðHTm ðxÞ½P 1mm  UHm ðtÞ þ 2tÞ ¼ HTm ðxÞ½P1mm  UP1=2 ¼ J 1=2 t : mm H m ðtÞ þ 1=2 @t Cð5=2Þ @t

ð47Þ

Substituting Eqs. (46) and (47) into Eq. (42), we have T

T

1 T 1 1=2 HTm ðxÞ½P2=3 mm  UP mm H m ðtÞ þ H m ðxÞ½P mm  UP mm H m ðtÞ ¼ 0:

ð48Þ

According to Eq. (48), we may find that U ¼ 0 is the solution of Eq. (48). Substituting U ¼ 0 into Eq. (45), we get uðx; tÞ ¼ x2 þ t2 which is the exact solution of the initial fractional partial differential equation. Example 3. Consider this equation

@au @bu þ ¼ f ðx; tÞ; @xa @tb

0 6 x; t 6 1;

ð49Þ

   ¼ 2x; @u such that @u ¼ 2t; uð0; tÞ ¼ t 2 þ 1; @x t¼0 @t x¼0 tion is ðx2 þ 1Þðt2 þ 1Þ. @2 u Let @x@t  HTm ðxÞUHm ðtÞ, then we have

@u ¼ @t

Z

t

0

Z 0

x

Cð3aÞ

þ

Cð3Þðx2 þ1Þt2b Cð3bÞ

, the exact solu-

  Z t @2u @u @u dt þ   ½HTm ðxÞUHm ðtÞdt þ  ¼ HTm ðxÞUP1mm Hm ðtÞ þ 2x; @x@t @x t¼0 @x t¼0 0

ð50Þ

  Z x T @2u @u @u dx þ   ½HTm ðxÞUHm ðtÞdx þ  ¼ HTm ðxÞ½P1mm  UHm ðtÞ þ 2t: @x@t @t x¼0 @t 0 x¼0

ð51Þ

1 0.8 0.6 u(x,t)

@u ¼ @x

uðx; 0Þ ¼ x2 þ 1, and f ðx; tÞ ¼

Cð3Þx2a ðt2 þ1Þ

0.4 0.2 0 1 0.5

t

0

0

0.2

0.4

x

Fig. 2. Numerical solution of m ¼ 16.

0.6

0.8

1

73

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

1

u(x,t)

0.8 0.6 0.4 0.2 0 1 0.5

t

0

0

0.2

0.4

0.6

0.8

1

x

Fig. 3. Numerical solution of m ¼ 32.

1

u(x,t)

0.8 0.6 0.4 0.2 0 1 0.5

t

0

0

0.2

0.4

0.6

0.8

1

x

Fig. 4. Exact solution for Example 1.

Therefore T

T

uðx; tÞ  HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ x2 þ uð0; tÞ ¼ HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ x2 þ t 2 þ 1:

ð52Þ

Then we can get

  h iT @au Cð3Þ 2a 1a @u a 1  J 1a ðHTm ðxÞUP1mm Hm ðtÞ þ 2xÞ ¼ HTm ðxÞ P1 ¼ J x ; mm UP mm Hm ðtÞ þ a @x @x Cð3  aÞ

ð53Þ

  T T @bu @u Cð3Þ 2b  J 1b ðHTm ðxÞ½P1mm  UHm ðtÞ þ 2tÞ ¼ HTm ðxÞ½P1mm  UP1b ¼ J 1b t : mm H m ðtÞ þ b @t Cð3  bÞ @t

ð54Þ

Substituting Eqs. (53) and (54) into Eq. (49), we have T

T

a 1 T 1 1b HTm ðxÞ½P1 mm  UP mm Hm ðtÞ þ H m ðxÞ½P mm  UP mm H m ðtÞ ¼ gðx; tÞ; Cð3Þx2a t2 Cð3aÞ

ð55Þ

Cð3Þx2 t 2b Cð3bÞ .

where gðx; tÞ ¼ þ Similarly, gðx; tÞ can be expressed as follows

gðx; tÞ  HTm ðxÞGHm ðtÞ;

ð56Þ

where G ¼ ½g ij mm . Dispersing Eqs. (55) and (56) by the points ðxi ; t j Þ, i ¼ 1; 2; . . . ; m and j ¼ 1; 2; . . . ; m, we can obtain T

T

a 1 T 1 1b T HT ½P1 mm  UP mm H þ H ½P mm  UP mm H ¼ H GH:

ð57Þ

74

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

Namely T

T

T

1 1a 1b 1 1 ½P 1 mm  ½P mm  U þ UP mm P mm ¼ ½P mm  GP mm :

ð58Þ

Eq. (58) is a Sylvester equation. We obtain U by solving Eq. (58). Then using Eq. (52), we get the numerical solution of uðx; tÞ. Taking a ¼ 1=2; b ¼ 1=3, we can achieve the absolute errors for different m. The absolute errors are shown in Table 1. From the Table 1, we can see clearly that the absolute errors become more and more small when m increases. The numerical results and the exact result for x ¼ 0:25, m ¼ 64 are shown in Fig. 5. From the Fig. 5, we find easily that the numerical solutions are in good agreement with the exact solution. Example 4. Consider the below fractional partial differential equation

@au @bu þ ¼ cos x þ cos t; @xa @tb

0 6 x; t 6 1;

   ¼ cos x; @u subject to @u ¼ cos t; @x t¼0 @t x¼0 2 @ u Let @x@t  HTm ðxÞUHm ðtÞ, then

@u ¼ @x @u ¼ @t

Z

t

0

Z

x

0

ð59Þ

uð0; tÞ ¼ sin t;

uðx; 0Þ ¼ sin x.

  Z th i @2u @u @u dt þ   HTm ðxÞUHm ðtÞ dt þ  ¼ HTm ðxÞUP1mm Hm ðtÞ þ cos x; @x@t @x t¼0 @x t¼0 0

ð60Þ

  Z x T @2u @u @u dx þ   ½HTm ðxÞUHm ðtÞdx þ  ¼ HTm ðxÞ½P1mm  UHm ðtÞ þ cos t: @x@t @t x¼0 @t x¼0 0

ð61Þ

Hence we have T

T

uðx; tÞ  HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ sin x þ uð0; tÞ ¼ HTm ðxÞ½P1mm  UP1mm Hm ðtÞ þ sin x þ sin t:

ð62Þ

Substituting Eqs. (60) and (61) into Eq. (59) when a ¼ b ¼ 1 T

HTm ðxÞUP1mm Hm ðtÞ þ HTm ðxÞ½P1mm  UHm ðtÞ ¼ 0:

ð63Þ

U ¼ 0 is the exact solution of Eq. (63). We can get uðx; tÞ  sin x þ sin t by using Eq. (62). When a ¼ b ¼ 1, the exact solution of initial partial differential equation is sin x þ sin t. When a; b–1, we have

 

 @au 1a @u a T 1 1a  J 1a HTm ðxÞUP1mm Hm ðtÞ þ cos x ¼ HTm ðxÞ½P 1 ¼ J ðcos xÞ; mm  UP mm H m ðtÞ þ J @xa @x

ð64Þ

  T T @bu 1b @u 1b  J 1b ðHTm ðxÞ½P1mm  UHm ðtÞ þ cos tÞ ¼ HTm ðxÞ½P1mm  UP1b ¼ J ðcos tÞ: mm H m ðtÞ þ J @t @tb

ð65Þ

Substituting Eqs. (64) and (65) into Eq. (59), we have T

T

a 1 T 1 1b HTm ðxÞ½P1 mm  UP mm H m ðtÞ þ H m ðxÞ½P mm  UP mm H m ðtÞ ¼ gðx; tÞ;

where gðx; tÞ ¼ cos x  J

Let cos x 

1a

ðcos xÞ þ cos t  J

uT1 Hm ðxÞ;

cos t 

1b

ð66Þ

ðcos tÞ.

uT2 Hm ðtÞ;

T

where u1 ¼ ½c0 ; c1 ; . . . ; cm1  ; u2 ¼ ½c00 ; c01 ; . . . ; c0m1 T . Then gðx; tÞ will be a T 1b gðx; tÞ ¼ cos x  uT1 P1 mm H m ðxÞ þ cos t  u1 P mm Hm ðtÞ:

ð67Þ

Table 1 The absolute errors of different m for Example 3. ðx; tÞ

m¼8

m ¼ 16

m ¼ 32

m ¼ 64

(0, 0) (1/8, 1/8) (2/8, 2/8) (3/8, 3/8) (4/8, 4/8) (5/8, 5/8) (6/8, 6/8) (7/8, 7/8)

5.415112e006 2.774770e005 5.600267e005 4.902546e005 7.292534e005 3.242363e005 6.928689e005 7.414284e005

3.424102e007 3.469165e006 4.710846e006 4.698928e006 4.491834e006 4.562334e006 4.998009e006 5.722911e006

2.164048e008 3.051360e007 3.149939e007 3.539285e007 4.675729e007 6.168849e007 7.817335e007 9.562386e007

1.366948e009 2.210144e008 3.298079e008 5.506236e008 8.024396e008 1.074419e007 1.363378e007 1.667197e007

75

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76 2.4

Numerical solution Exact solution

2.2

u(0.25,t)

2 1.8 1.6 1.4 1.2 1

0

0.1

0.2

0.3

0.4

0.5

t

0.6

0.7

0.8

0.9

1

Fig. 5. Comparison of numerical solutions and exact solution x ¼ 0:25; m ¼ 64.

2

u(x,t)

1.5 1 0.5 0 1 0.5

t

0

0

0.2

0.4

0.6

0.8

1

x

Fig. 6. Numerical solution of a ¼ 3=4; b ¼ 2=3.

2

u(x,t)

1.5 1 0.5 0 1 0.5

t

0

0

0.2

0.4

0.6

x

Fig. 7. Numerical solution of a ¼ 3=5; b ¼ 1=3.

0.8

1

76

L. Wang et al. / Applied Mathematics and Computation 227 (2014) 66–76

Similarly, gðx; tÞ can be also expressed as follows

gðx; tÞ  HTm ðxÞGHm ðtÞ;

ð68Þ

where G ¼ ½g ij mm . Dispersing Eqs. (66) and (68) by the points ðxi ; tj Þ, i ¼ 1; 2; . . . ; m and j ¼ 1; 2; . . . ; m, we have T

T

a 1 T 1 1b T HT ½P 1 mm  UP mm H þ H ½P mm  UP mm H ¼ H GH:

ð69Þ

Thus T

T

T

1 1a 1b 1 1 ½P 1 mm  ½P mm  U þ UP mm P mm ¼ ½P mm  GP mm :

ð70Þ

Eq. (70) is a Sylvester equation. We can get U by solving Eq. (70). Then applying Eq. (62), we obtain the approximation of uðx; tÞ. Figs. 6 and 7 show the numerical solutions for different values of a; b. Here, we may take m ¼ 32. Compared with the generalized differential transform method in Ref. [19], taking advantage of above technique greatly reduces computation. What is more, the method in this paper is easy implementation.

5. Conclusion The aim of this paper is to develop an effective and accurate method for solving fractional partial differential equations. We give the Haar wavelet operational matrix of fractional order integration. This matrix is used to reduce the fractional partial differential equations to Sylvester equations. The convergence analysis of the Haar wavelet bases is given in Section 3. The numerical solutions obtained using the suggested method show that numerical solutions are in very good coincidence with the exact solution. Acknowledgment The International Science &Technology Cooperation Program of China (No. 2012DFG61930). References [1] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and applications of the fractional differential equations, Mathematical Studies, vol. 204, Elsevier, Amsterdam, 2006. [2] K.B. Oldham, Fractional differential equations in electrochemistry, Advances in Engineering Software 41 (2010) 9–12. [3] R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific Publishing Company, Singapore, 2000. [4] H.H. Sun, A.A. Abdelwahad, B. Onaral, Linear approximation of transfer function with a pole of fractional order, IEEE Transactions on Automatic Control 29 (1984) 441–444. [5] M. Ichise, Y. Nagayanagi, T. Kojima, An analog simulation of noninteger order transfer functions for analysis of electrode process, Journal if Electroanalysis Chemistry 33 (1971) 253–265. [6] R.C. Koeller, Applications of fractional calculus to the theory of viscoelasticity, Journal of Applied Mechanics 51 (1984) 299–307. [7] I.L. EI-Kalla, Convergence of the Adomian method applied to a class of nonlinear integral equations, Applied Mathematics and Computation 21 (2008) 372–376. [8] M.M. Hosseini, Adomian decomposition method for solution of nonlinear differential algebraic equations, Applied Mathematics and Computation 181 (2006) 1737–1744. [9] S. Momani, Z. Odibat, Generalized differential transform method for solving a space and time-fractional diffusion-wave equation, Physics Letters A 370 (2007) 379–387. [10] Z. Odibat, S. Momani, Generalized differential transform method: application to differential equations of fractional order, Applied Mathematics and Computation 197 (2008) 467–477. [11] Z. Odibat, A study on the convergence of variational iteration method, Mathematical and Computer Modelling 51 (2010) 1181–1192. [12] H. Sun, W. Chen, C. Li, et al, Finite difference schemes for variable-order time fractional diffusion equation, International Journal of Bifurcation and Chaos 22 (4) (2012) 1250085. [13] M.M. Meerschaert, H.P. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, Journal of Computational Physics 211 (2006) 249–261. [14] I. Hashim, O. Abdulaziz, Homotopy analysis method for fractional IVPs, Communications in Nonlinear Science and Numerical Simulation 14 (3) (2009) 674–684. [15] Y.M. Chen, Y.B. Wu, et al, Wavelet method for a class of fractional convection–diffusion equation with variable coefficients, Journal of Computational Science 1 (3) (2010) 146–149. [16] H. Jafari, S.A. Yousefi, Application of Legendre wavelets for solving fractional differential equations, Computers and Mathematics with Application 62 (3) (2011) 1038–1045. [17] Y.M. Chen, M.X. Yi, C.X. Yu, Error analysis for numerical solution of fractional differential equation by Haar wavelets method, Journal of Computational Science 5 (3) (2012) 367–373. [18] I. Podlubny, Fractional Differential Equations, Academic press, 1999. [19] Zaid Odibat, Shaher Momani, A generalized differential transform method for linear partial differential equations of fractional order, Applied Mathematics Letters 21 (2008) 194–199. [20] Y. Zhang, A finite difference method for fractional partial differential equation, Applied Mathematics and Computation 215 (2009) 524–529. [21] C.F. Chen, C.H. Hsiao, Haar wavelet method for solving lumped and distributed-parameter systems, IEE Proceedings Control Theory and Applications 144 (1) (1997) 87–94. [22] Y.L. Li, Haar wavelet operational matrix of fractional order integration and its applications in solving the fractional order differential equations, Applied Mathematics and Computation 216 (8) (2010) 2276–2285.