Obtaining artificial boundary conditions for fractional sub-diffusion equation on space two-dimensional unbounded domains

Computers and Mathematics with Applications 68 (2014) 13–26

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Obtaining artificial boundary conditions for fractional sub-diffusion equation on space two-dimensional unbounded domains Rezvan Ghaffari, S. Mohammad Hosseini ∗ Department of Applied Mathematics, Faculty of Mathematical Sciences, Tarbiat Modares University, P.O. Box 14115-175, Tehran, Iran

article

info

Article history: Received 11 November 2013 Received in revised form 24 March 2014 Accepted 3 May 2014 Available online 21 May 2014 Keywords: Fractional sub-diffusion equation Artificial boundary condition Two-dimensional unbounded domain Laplace transform Modified Bessel equation

abstract In this paper, we consider fractional sub-diffusion equation on two-dimensional unbounded domains and obtain an exact and some approximating artificial boundary conditions for it by a joint application of the Laplace transform and the Fourier series expansion. In order to test the derived artificial boundary condition, after reducing the main problem to an initial–boundary value problem on bounded domain by imposing the obtained boundary condition, we just use classical Crank–Nicolson method for space variables and L1 approximation for the fractional time derivative. Some numerical examples are given which confirm the effectiveness of the proposed artificial boundary conditions. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction Fractional differential equations and fractional calculus arise in various application problems in physical, biological, geological, and chemical systems; mechanical engineering, environment sciences, signal processing, systems identification, electrical and control theory, etc., see [1–6]. The fractional sub-diffusion equation is a subclass of anomalous diffusive systems, which can be obtained from the standard parabolic PDEs by replacing the first-order time derivative with a fractional derivative of order α, 0 < α < 1. For some general results regarding the existence and uniqueness of the solution of fractional differential equations we can refer the reader to [1,7]. We consider fractional sub-diffusion equation on two-dimensional unbounded domains and obtain an exact and some approximating artificial boundary conditions for it by a joint application of the Laplace transform and the Fourier series expansion. As far as we searched the related literature there were no reports on similar results as we are investigating in this paper. The problem of fractional sub-diffusion equation on bounded domains in two space dimensions has been widely studied from different aspects. There are many numerical methods reported for related problems on bounded domain. For example, explicit and implicit finite difference schemes [8,9], implicit RBF meshless approach [10], alternating direction implicit (ADI) scheme [11,12], high-order compact finite difference [13], orthogonal spline collocation (OSC) [14], Kansa-type method [15], etc. However, considering such problems on unbounded two-dimensional domains could be very interesting and important in applications. The great difficulties to obtain the numerical solutions of such problems on unbounded physical domains lie

∗

Corresponding author. E-mail addresses: [email protected] (R. Ghaffari), [email protected], [email protected] (S.M. Hosseini).

http://dx.doi.org/10.1016/j.camwa.2014.05.006 0898-1221/© 2014 Elsevier Ltd. All rights reserved.

14

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

in the unboundedness of physical domains. There are a lot of papers that have discussed the analytical solutions of fractional problems. For some sort of fractional equations on bounded domains, for example see [16–18] and on unbounded domains using various methods such as Laplace and Fourier transforms, for example see [19–21]. Solutions of such problems are often represented in terms of Fox-H functions, M-Wright functions, Green functions and Mittag-Leffler functions, where these functions are in the form of series and their computation is not so easy. In this work we focus on a novel derivation of artificial boundary conditions, because the fractional problems on unbounded domains in application usually need to be solved numerically with efficient use of any suitable numerical method. However, the concept of uniqueness and continuous dependence of the solution of the fractional sub-diffusion equation with a polar coordinate in two-dimensional unbounded domains, will be investigated in detail in future works. Artificial boundary method (ABM) has been one of the most important and efficient methods for numerical solution of partial differential equations (PDEs) in unbounded domains and has been widely applied for integer order equations on unbounded domains in problems such as wave equations [22–24], parabolic equations [25–27], more different types of equations [28] and one-dimensional fractional sub-diffusion equations [29–31]. In the artificial boundary method, at first, one divides the original unbounded domain by an artificial boundary into subdomains: a bounded domain and an unbounded residual domain. Then the problem in unbounded residual domain is considered for deriving the exact and approximating artificial boundary conditions on the artificial boundary. Second, the obtained artificial boundary condition is used by any suitable numerical methods, newly introduced or from those already reported in the literature over bounded domains, to solve the problem in bounded domain. So, in this paper, we investigate the derivation of artificial boundary condition for the fractional sub-diffusion equation on space two-dimensional unbounded domains. The results should be very interesting for applications in sciences and engineering. Using the obtained artificial boundary conditions, any of the already introduced numerical methods by various authors for bounded two-dimensional domains can easily be applied to the problem on bounded subdomain. However, to test our findings we just apply the Crank–Nicolson finite difference method without losing net worth of the obtained results. The discussion in the rest of paper will be as follows. In Section 2, we express the problem in unbounded domain and introduce an artificial boundary. In Section 3, we derive artificial boundary conditions by using the Laplace transformation and some interesting properties of the solutions of the modified Bessel equations. Finally, in Section 4, some numerical results are presented to show the effectiveness of the results proposed in this paper. 2. The problem Let D ⊂ R2 denote a bounded domain where D = x ∈ R2 | ∥x∥ < a . Suppose



Γ0T = ∂ D × [0, T ],

¯, Dc = R2 \D



ΩcT = Dc × [0, T ].

Consider the following initial–boundary value problem: c α 0 Dt u

(x, t ) − 1u(x, t ) = f (x, t ),

(x, t ) ∈ ΩcT ,

u|Γ0 = g (x, t ),

(x, t ) ∈ Γ0T ,

(2)

u|t =0 = u0 (x),

x∈D ,

(3)

u(x, t ) → 0,

c

when ∥x∥ → +∞, t > 0,

(4)

c α 0 Dt

(0 < α < 1) is the Caputo fractional derivative of order α defined by [1]  t ∂ u(x, λ) 1 c α (t − λ)−α dλ, D u ( x , t ) = 0 t Γ (1 − α) 0 ∂λ  and x = (x1 , x2 ) is the Cartesian coordinate in space, ∥x∥ = x21 + x22 . ∆ is the Laplacian operator, i.e., ∆ =

where

∂2

(1)

∂2

∂2 ∂ x21

+

∂2 ∂ x22

=

+ 1r ∂∂r + r12 ∂θ 2 . Functions f (x, t ), g (x, t ) and u0 (x) are given smooth functions, and f (x, t ) and u0 (x) vanish outside a ∂r2 disc with the radius b, namely, f (x, t ) = 0,

u0 (x) = 0,

for ∥x∥ > b. We introduce an artificial boundary

  Γ = (x, t )|x ∈ Dc , ∥x∥ = b, 0 6 t 6 T ,

(5)

where a < b. Γ divides the domain ΩcT into two parts,

  ΩiT = (x, t )|x ∈ Dc , ∥x∥ < b, 0 6 t 6 T , ΩeT = {(x, t )| ∥x∥ > b, 0 6 t 6 T } .

(6)

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

15

3. Derivation of the artificial boundary conditions The aim of paper is to derive the exact and some approximating artificial boundary conditions. First, we consider the problem (1)–(4) on the unbounded subdomain ΩeT in the polar coordinate, that is, u(r , θ , t ) satisfies

∂ 2u 1 ∂ u 1 ∂ 2u + + , (r , θ , t ) ∈ ΩeT , (7) ∂r2 r ∂r r 2 ∂θ 2 u|Γ = u(b, θ , t ), (8) u| t = 0 = 0 , (9) u(r , θ , t ) → 0, when r → +∞, t > 0, (10) where ΩeT = {r > b, θ ∈ [0, 2π ], t ∈ [0, T ]}. Since u(b, θ , t ) is an unknown function, problem (7)–(10) is an incompletely posed problem. It cannot be solved independently. By assuming that the boundary function u(b, θ , t ) is given, we try to find the solution of problem (7)–(10). Let us express the function u(b, θ , t ) by Fourier series: c α 0 Dt u

(r , θ , t ) =

u(b, θ , t ) = a0 (t ) +

∞ 

(an (t ) cos nθ + bn (t ) sin nθ ) ,

(11)

n =1

a0 (t ) = an (t ) = bn (t ) =

2π 1

π 1

2π



1

u(b, θ , t )dθ ,

(12)

0 2π

 0



u(b, θ , t ) cos nθ dθ , (13)

2π

u(b, θ , t ) sin nθ dθ , π 0 where n = 1, 2, . . . . Similarly, the Fourier series of the solution function u(r , θ , t ) can be written as u(r , θ , t ) = u0 (r , t ) +

∞ 

{un (r , t ) cos nθ + vn (r , t ) sin nθ} ,

(14)

n=1

where its Fourier coefficients can also be written similar as above. By substituting (14) into (7), we obtain

  ∞   1 ∂ u0 1 ∂ un n2 ∂ 2 u0 ∂ 2 un c α − + − + 2 un cos nθ (r , t ) − 0 Dt un (r , t ) − ∂r2 r ∂r ∂r2 r ∂r r n=1    1 ∂vn n2 ∂ 2 vn + c0 Dαt vn (r , t ) − − + 2 vn sin nθ = 0. 2 ∂r r ∂r r

c α 0 Dt u0

The above equation holds if and only if the coefficients of the series vanish, hence we have: (I) Functions un (r , t ), n = 0, 1, 2, . . . , satisfy 1 ∂ un n2 ∂ 2 un + un , r > b , − ∂r2 r ∂r r2 un |r =b = an (t ), un |t =0 = 0, un (r , t ) → 0, when r → +∞, t > 0. (II) Functions vn (r , t ), n = 1, 2, . . . , satisfy c α 0 Dt un

(r , t ) =

(15) (16) (17) (18)

∂ 2 vn 1 ∂vn n2 + − 2 vn , r > b, (19) 2 ∂r r ∂r r vn |r =b = bn (t ), (20) vn |t =0 = 0, (21) vn (r , t ) → 0, when r → +∞, t > 0. (22) Then, if one obtains un (r , t ), n = 0, 1, 2, . . . , and vn (r , t ), n = 1, 2, . . . , respectively, as the solutions of the problems (15)–(18) and (19)–(22), the function u(r , θ , t ), as the solution of problem (7)–(10) given by (14), can be obtained consec α 0 Dt n

v (r , t ) =

quently. Now we start to solve those problems by applying the Laplace transformation, the details of which are given here only for problem (15)–(18) as it applies similarly for other problem. Let us denote the Laplace transform of function un as u˜ n , which

16

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

is defined as u˜ n (r , s) =

+∞



un (r , λ)e−sλ dλ.

0

The Laplace transform of (15)–(18) and the Caputo fractional derivatives, see [1,32], such as

L[c0 Dαt f (t )](s) = sα L[f (t )](s) − sα−1 f (0),

0 < α 6 1,

lead to sα u˜ n (r , s) =

d2 u˜ n (r , s) dr 2

+

1 du˜ n (r , s) r

dr

−

n2 r2

u˜ n (r , s),

r > b,

(23)

u˜ n (r , s)|r =b = a˜ n (s),

(24)

u˜ n (r , s) → 0,

(25)

when r → +∞.

Then Eqs. (23)–(25) can be written as d2 u˜ n (r , s) dr 2

+

1 du˜ n (r , s) r

dr 1

u˜ n (r , s)|r =b =

r

r > b,

(26)

· sa˜ n (s),

s

u˜ n (r , s) → 0,

  n2 − sα + 2 u˜ n (r , s) = 0,

(27)

when r → +∞.

(28)

Obviously, the solution of Eq. (26) can then be easily expressed in the form

˜ n (r , s), u˜ n (r , s) = [sa˜ n (s)]G

(29)

˜ n (r , s) is the solution of boundary value problem where G ˜ n ( r , s) d2 G dr 2

+

˜ n (r , s) 1 dG r

dr 1

˜ n (r , s)|r =b = G

s

˜ n (r , s) → 0, G α



α

− s +

n2 r2



˜ n (r , s) = 0, G

r > b,

(30)

,

(31)

when r → +∞.

(32)

Eq. (30), with s = β , changes to

˜ n ( r , s) d2 G dr 2

2

+

˜ n (r , s) 1 dG r

˜ n (r , s)|r =b = G

dr 1 s

˜ n (r , s) → 0, G

  n2 − β 2 + 2 G˜ n (r , s) = 0, r

r > b,

(33)

,

(34)

when r → +∞.

(35)

Eq. (33) is known as the modified Bessel equation of order n [33–38]. The equation has two independent solutions In (β r ) and Kn (β r ), respectively, known as the first and second kind modified Bessel functions. Hence the general solution can be written as

˜ n (r , s) = A(β)In (β r ) + B(β)Kn (β r ). G

(36)

From the infinity condition (35), we must have A(β) = 0, because from [33–38]

 (µ − 1)(µ − 9) (µ − 1)(µ − 9)(µ − 25) 1− − + ··· , In (z ) ∼ √ + 8z 2!(8z )2 3!(8z )3 2π z    π −z µ − 1 (µ − 1)(µ − 9) (µ − 1)(µ − 9)(µ − 25) Kn (z ) ∼ e 1+ + + + ··· , 2z 8z 2!(8z )2 3!(8z )3 ez



µ−1

for |z | large and µ = 4n2 . Satisfying the boundary condition (34) yields B(β) =

1 sKn (β b)

.

|arg z | <

π 2

|arg z | <

, 3π 2

,

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

17

Therefore, the solution of problem (33)–(35) is as follows: Kn (β r )

˜ n (r , s) = G

sKn (β b)

,

r > b.

(37)

˜ n (r , s), see [39,40], is given by The function Gn (r , t ), as the Laplace inverse transform of G Gn (r , t ) =

c +i∞



1 2π i

˜ n (r , s)ds, est G

(38)

c −i ∞

˜ n (r , s) lie to the left of this line. Subject to the integration being along a line from c − i∞ to c + i∞ such that all poles of G the this requirement the value of c is arbitrary. The integral (38) will be solved by the method of residues. ˜ n (r , s). The modified Bessel function of the second kind can be Let us first recall some facts about the function G expressed as n −1 1

Kn (z ) =

2 r =0

(−1)r

  ∞  (n − r − 1)! (z /2)n+2r 1 (2/z )n−2r + (−1)n+1 ln(z /2) − [ψ(r + 1) + ψ(n + r + 1)] , r! r !(n + r )! 2 r =0

where

ψ(n + r + 1) = 1 + 1/2 + 1/3 + · · · + 1/(n + r ) − γ , ˜ n (r , s) has a branch point at the origin. and γ = 0.5772156649015329 . . . . So, due to the logarithmic term, Kn (z ) and hence G On the other hand, since according to [33,41–43],all the finite zeros of Kn (z ) for n ≥ 2 are complex numbers with negative real parts Re(z ) < 0 and cannot lie on the negative real axis, the branch cut is considered to lie in the quadrants for which π < |arg z | < π . So, in other words, G˜ n (r , s) has poles for n ≥ 2 which are known to occur in complex conjugate pairs and 2 has no pole for n = 0, 1, as Kn (z ) has no finite zeros for such values of n. The complex zeros of Kn (z ) are given in [41,42]. Now to calculate the integral, see [39,40], we take the closed contour consisting of: z = Reiθ , π /2 < θ < π ; arg(z ) = π , ϵ < |z | < R; z = ϵ eiθ , −π < θ < π ; arg(z ) = −π , ϵ < |z | < R; z = Reiθ , −π < θ < −π /2; z = ξ − iR, 0 < ξ < c; z = c + iξ , − R < ξ < R; z = ξ + iR, 0 < ξ < c, with R large enough to include all the poles and ϵ sufficiently small. It can be shown that by the limiting property of the integrand, the integrals over those pieces of the contour which depend on R vanish in the limits as R −→ ∞. Hence (38) turns by Cauchy’s Theorem to 

Gn (r , t ) =



1



˜ n (r , s); ak + Res est G

k



2π i

˜ n (r , s)ds, est G L

where ‘‘Res’’ denotes the residue and



˜ n (r , s)ds = e G st

L

∞



e

−st



ε

π

  ˜Gn (r , se−iπ ) − G˜ n (r , s+iπ ) ds + i

iθ

˜ n (r , εeiθ )dθ , et εe (ε eiθ )G

−π

where ε is any sufficiently small positive number. Hence, we can write [40]



Gn (r , t ) =

k

 ∞    1 ˜ ˜ n (r , se−iπ ) − G˜ n (r , se+iπ ) ds + An , Res e Gn (r , s); ak + e−st G 2π i 0 

st

(39)

˜ n (r , s). where An = lims→0 sG The first term on the right-hand side of (39), based on the particular distribution of the poles, is zero, (see Theorems 1 and 2 in pp. 284 and 285 in [44]). The second term on the right-hand side of (39) can be expressed as: 1 2π i

∞



−st

e



˜ n (r , se G

−iπ

0

 ∞ −st  1 e ) − G˜ n (r , se ) ds = 2π i 0 s =

For the third term in (39), we obtain

˜ n (r , s) = lim An = lim sG s→0

β→0

Kn (β r ) Kn (β b)

.

From [33–38], while z → 0 Kn (z ) ∼



iπ

 −Ln(z ), 2n−1 (n − 1)!z −n ,

n = 0, n = 1, 2, . . . .

1

απ i

∞

 0

e−β

β

α



α

−

Kn β re 2 π i



Kn β be 2 π i



2

αt









α

−

Kn β be 2 π i





α

 ds

Kn β be− 2 π i

α

Kn β re 2 π i

α

Kn β re− 2 π i



α



α

 dβ.

Kn β re− 2 π i



Kn β be− 2 π i



(40)

18

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

Hence A0 = lim

β→0

An = lim

β→0

K0 (β r ) K0 (β b) Kn (β r ) Kn (β b)

= lim

β→0

= lim

β→0

Ln(β r ) Ln(β b)

= 1,

2n−1 (n − 1)!(β r )−n 2n−1 (n − 1)!(β b)−n

 n b

=

r

.

So, we have

 n An =

b r

,

n = 0, 1, 2, . . . .

(41)

Therefore, from (40), (41), Eq. (39) finally becomes

 n b

Gn (r , t ) =

+

r

1

απ i

∞

 0

e−β

2

αt



β

α



α

−

Kn β re 2 π i



Kn β be 2 π i



α



α

 dβ.

Kn β re− 2 π i



Kn β be− 2 π i



(42)

By applying convolution theorem of the Laplace transform and by using (29) and L−1 sa˜ n (s) =



un (r , t ) = L−1



  da (t ) n ∗ Gn (r , t ) = (sa˜ n (s)) · G˜ n (r , s) = dt

t 0

dan (λ) dλ

Gn (r , t − λ)dλ,



dan (t ) , dt

we then have (43)

where un (r , t ) is the solution of problem (15)–(18). Furthermore, we have

   t dan (λ) ∂ Gn (r , t − λ)  ∂ un (r , t )   dλ. =  ∂ r r =b dλ ∂r 0 r =b

(44)

On the other hand, from (42),

  − α πi   ∞ −β α2 t ∂ Kn (β re α2 π i ) ∂ Kn (β re 2 ) | r =b | r =b −n ∂ Gn (r , t )  1 e ∂ r ∂ r   dβ. = + − α α  ∂r b απ i 0 β Kn (β be 2 π i ) Kn (β be− 2 π i ) r =b

(45)

From [33–38]

 −K (z ), ∂ Kn (z )  1 = −1  ∂z [Kn−1 (z ) + Kn+1 (z )] , 2

n = 0, n = 1, 2, . . . .

(46)

Hence, when n = 0

     α π i  ∞ − α πi 2 1 ∂ G0 (r , t )  α π i K1 β be 2 α π i K1 β be 2 −β α t − 2 = e dβ, −e 2   α  +e α   ∂r απ i 0 K0 β be 2 π i K0 β be− 2 π i r =b

(47)

and, when n = 1, 2, . . .

    α π i α   ∞ 2 + Kn+1 β be 2 π i ∂ Gn (r , t )  −n 1 α π i Kn−1 β be 2 α −β t 2 + = e −e  α   ∂r b 2απ i 0 Kn β be 2 π i r =b    α  − α2 π i + Kn+1 β be− 2 π i − α2 π i Kn−1 β be +e dβ.  α  Kn β be− 2 π i

(48)

In other words, (47)–(48) can be expressed as:

 ∂ Gn (r , t )   ∂r

r =b

=

−n b

+ Rn (t ),

where

     α π i  ∞ − α πi 2 1 α π i K1 β be 2 α π i K1 β be 2  α −β t −  2 e −e 2 dβ, n = 0,    α  +e α    απ i 0 K0 β be 2 π i K0 β be− 2 π i        α  α   1  ∞ 2 Kn−1 β be 2 π i + Kn+1 β be 2 π i α −β α t Rn (t ) = e −e 2 π i   α  2απ i 0 Kn β be 2 π i        α   − α πi  + Kn+1 β be− 2 π i α π i Kn−1 β be 2  −  +e 2 dβ, n = 1, 2, . . . .   α  Kn β be− 2 π i

(49)

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

19

Combining (44) and (49), yields the following result for the derivative of function un (r , t ) at r = b:

   t ∂ un (r , t )  dan (λ) ∂ Gn (r , t − λ)   dλ =  ∂ r  r =b dλ ∂r 0 r =b  t  t dan (λ) dan (λ) −n dλ + Rn (t − λ)dλ = dλ b dλ 0 0  t ∂ un (b, λ) −n un (b, t ) + Rn (t − λ)dλ. = b ∂λ 0

(50)

Then, by applying a similar procedure to the problem (19)–(22) the following results are also obtained for function vn (r , t ):

vn (r , t ) =

t



dbn (λ)

Gn (r , t − λ)dλ,

(51)

  t ∂vn (r , t )  ∂vn (b, λ) −n v ( b , t ) + Rn (t − λ)dλ. = n ∂ r r =b b ∂λ 0

(52)

dλ

0

and

Now we are ready to use the obtained results (50) and (52) in Eq. (14) in which u0 (r , t ) = un (r , t ) =

vn (r , t ) =



1 2π 1 1

0 2π



π

2π

u(r , φ, t )dφ,

u(r , φ, t ) cos nφ dφ,

0



2π

u(r , φ, t ) sin nφ dφ π 0 where n = 1, 2, . . . . So, taking the derivative with respect to r from both sides of (14) and replacing (50) and (52), yields      ∞   ∂ u(r , θ , t )  ∂ u0 (r , t )  ∂ un (r , t )  ∂vn (r , t )  = + cos n θ + sin n θ  ∂r ∂r  ∂r  ∂r  r =b

r =b

=−

∞  n n =1

+

b

n =1

r =b

r =b

(un (b, t ) cos nθ + vn (b, t ) sin nθ )

 ∞  t  ∂ un (b, λ) ∂vn (b, λ) cos nθ + sin nθ Rn (t − λ)dλ. ∂λ ∂λ n =0 0

(53)

Finally, the exact boundary condition on the artificial boundary Γ is obtained as

  t  2π 1 ∂ u(b, φ, λ) ∂ u(r , θ , t )  = dφ R0 (t − λ)dλ  ∂r 2π 0 0 ∂λ r =b ∞   1  t 2π ∂ u(b, φ, λ) + cos n(φ − θ )dφ Rn (t − λ)dλ π n =1 0 0 ∂λ  2π ∞ 1  − n u(b, φ, t ) cos n(φ − θ )dφ. bπ n=1 0

(54)

In general, by considering the first few terms of the above summations, one obtains a series of approximating artificial boundary conditions on Γ such as 2π

∂ u(b, φ, λ) dφ R0 (t − λ)dλ ∂λ 0 0 r =b  N  1  t 2π ∂ u(b, φ, λ) + cos n(φ − θ )dφ Rn (t − λ)dλ π n =1 0 0 ∂λ  2π N 1  − n u(b, φ, t ) cos n(φ − θ )dφ bπ n=1 0   ∂u , = κ N u, ∂t where N = 0, 1, 2, . . . .   ∂ u(r , θ , t )   ∂r

=

1

 t

2π

(55)

20

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

Remark 3.1. For particular value α = 1, Eqs. (47) and (48) can be simplified as follows: Based on the details given in [33–38,45], we have

   1 1 −1  Kn ze 2 π i = π ie 2 nπ i [−Jn (z ) + iYn (z )], 2    Kn ze− 12 π i = − 1 π ie 12 nπ i [−Jn (z ) − iYn (z )], 2

and we obtain 1

Kn β re 2 π i

1

Kn β re− 2 π i

Jn (β r ) + iYn (β r ) −Jn (β r ) + iYn (β r ) −  −   = 1 1 π i − π i −Jn (β b) + iYn (β b) Jn (β b) + iYn (β b) Kn β be 2 Kn β be 2 −Jn (β r )Jn (β b) − iJn (β r )Yn (β b) + iYn (β r )Jn (β b) − Yn (β r )Yn (β b) = −J 2 n (β b) − Y 2 n (β b) Jn (β r )Jn (β b) + iYn (β r )Jn (β b) − iJn (β r )Yn (β b) + Yn (β r )Yn (β b) + −J 2 n (β b) − Y 2 n (β b)   Jn (β r )Yn (β b) − Jn (β b)Yn (β r ) . = 2i J 2 n (β b) + Y 2 n (β b)









Hence, we have Gn (r , t ) =

 n b

+

r

∞



2

π

e−β

2t



Jn (β r )Yn (β b) − Jn (β b)Yn (β r ) J 2 n (β b) + Y 2 n (β b)

β

0



dβ,

(56)

and by using the Wronskian relation W (Jn (z ), Yn (z )) = Jn (z )Yn′ (z ) − Jn′ (z )Yn (z ) =

2

πz

,

we obtain

  ∞ 2 ∂ Gn (r , t )  −n 4 e−β t   dβ, = − 2  ∂r b π b 0 β J 2 n (β b) + Y 2 n (β b) r =b that is the same result as that of [26] which discusses a similar problem only for the case α = 1.

(57) 

Remark 3.2. If we write (23)–(25) as the following: d2 u˜ n (r , s) dr 2

+

1 du˜ n (r , s) r

dr

  n2 − sα + 2 u˜ n (r , s) = 0, r

r > b,

(58)

u˜ n (r , s)|r =b = a˜ n (s),

(59)

u˜ n (r , s) → 0,

(60)

when r → +∞,

then, we can obtain another form of artificial boundary condition such as:

  t  2π ∂ u(r , θ , t )  1 = u(b, φ, λ)dφ R0 (t − λ)dλ  ∂r 2π 0 0 r =b ∞   1  t 2π + u(b, φ, λ) cos n(φ − θ )dφ Rn (t − λ)dλ, π n=1 0 0 where R0 (t ) =

1

∞



e

απ i



2

−β α t

−e

  α  α π i  β be− 2 π i α π i K1 β be 2 2 β α dβ,   + e2  α π i − α2 π i K0 β be K0 β be 2

− α2 π i K1

0

and Rn (t ) =

∞



1 2απ i α

+ e 2 πi

2

e



−β α t

−e

− α2 π i Kn−1

0 α



 α  α  β be− 2 π i + Kn+1 β be− 2 π i  −α  Kn β be 2 π i α π i 

Kn−1 β be 2 π i + Kn+1 β be 2







α π i

Kn β be 2



2

β α dβ. 

(61)

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

21

4. Numerical examples In order to demonstrate the effectiveness of the artificial boundary conditions given in this paper, some numerical examples are considered. We assume Γ0 = {x| ∥x∥ = a} where a < b, hence

ΩcT = {r > a, θ ∈ [0, 2π], t ∈ [0, T ]}, ΩiT = {a < r < b, θ ∈ [0, 2π ], t ∈ [0, T ]} , ΩeT = {r > b, θ ∈ [0, 2π ], t ∈ [0, T ]} . By rewriting the original problem (1)–(4) in polar coordinate c α 0 Dt u

(r , θ , t ) =

∂ 2u 1 ∂ u 1 ∂ 2u + + 2 2 + f (r , θ , t ), 2 ∂r r ∂r r ∂θ

(r , θ , t ) ∈ ΩcT ,

(62)

u|Γ0 = g (θ , t ),

(63)

u|t =0 = u(r , θ , 0),

(64)

u → 0,

(65)

when r → ∞,

we then use the obtained boundary conditions (55) to reduce the problem (62)–(65) to the following approximating problem on bounded domain ΩiT : c α 0 Dt u

(r , θ , t ) =

∂ 2u 1 ∂ u 1 ∂ 2u + + + f (r , θ , t ), ∂r2 r ∂r r 2 ∂θ 2

(r , θ , t ) ∈ ΩiT ,

(66)

u|Γ0 = g (θ , t ),

(67)

u|t =0 = u(r , θ , 0),

 ∂ u  ∂r 



= κ N u,

r =b

(68)

∂u . ∂t 

(69)

For the numerical solution of the problem (66)–(69) we use the classical Crank–Nicolson method for space variables and L1 approximation for the fractional time derivative to solve the related test examples, however, any other numerical methods recently used on bounded domain can be similarly applied. The intervals [0, T ], [a, b] and [0, 2π ] are divided into K , I and J equal parts, respectively, 0 = t0 < t1 < · · · < tK = T ,

a = r0 < r1 < · · · < rI = b,

0 = θ0 < θ1 < · · · < θJ = 2π .

Let I = J /6, τ = T /K , h = (b − a)/I, δ = 2π /J, and τ = h. For the temporal approximation, we use

 c α 0 Dt u

(ri , θj , tk ) = µ

ukij

−

k−1 

 (ak−n−1 −

)

ak−n unij

−

ak−1 u0ij

+ O(τ 2−α ),

(70)

n =1

as the discrete fractional derivative operator or L1 approximation [9], where µ = Γ (τ2−α) and an = (n + 1)1−α − n1−α . Now we have the following formula by using the classical Crank–Nicolson method for space variables and L1 approximation for the fractional time derivative −α

1 k−1 uki+1,j − 2uki,j + uki−1,j + uki+−11,j − 2uki,− j + ui−1,j

2h2

+

uki,j+1

−

2uki,j

+

uki,j−1

+

1 uki,− j +1

−

1 2uki,− j

2ri 2 δ 2

 −µ

uki,j

+

+

uki+1,j − uki−1,j + uki+−11,j − uki−−11,j

1 uki,− j−1

4hri k− 12

+ fi,j

 k−1  l 0 − (ak−l−1 − ak−l )ui,j − ak−1 ui,j = 0,

1 6 i 6 I, 1 6 j 6 J, 1 6 k 6 K ,

(71)

l =1

u0i,j = u(ri , θj , 0), uk0,j = g (θj , tk ), ukI+1,j − ukI−1,j 2h

1 6 i 6 I, 1 6 j 6 J,

(72)

1 6 j 6 J, 1 6 k 6 K ,

(73)

  ∂u = κN u(rI , ., .), (rI , ., .) (θj , tk ), ∂t

1 6 j 6 J, 1 6 k 6 K

(74)

22

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26 k−1/2

where uki,j = u(ri , θj , tk ), fi,j these integrals: tk



∂ u(b, φ, λ) cos n(φ − θj )dφ Rn (tk − λ)dλ, ∂λ

0

0 2π



2π



= 12 (fi,kj−1 + fi,kj ) and uki,0 = uki,J . Due to the operator κN , at each time step we need to compute

u(b, φ, λ) cos n(φ − θj )dφ,

n = 0, 1, 2, . . . , N ,

(75)

n = 1, 2, . . . , N .

(76)

0

It is easy to calculate the integrals in (76), 2π



u(b, φ, λ) cos n(φ − θj )dφ

0

=

J −1   s=0

=

θs+1

θs

J −1 1 

n2

δ

 ukI,s

+

(ukI,s+1 − ukI,s )(φ − θs )



δ

cos n(φ − θj )dφ

ukI,s 2 cos n(θs − θj ) − cos n(θs+1 − θj ) − cos n(θs−1 − θj ) .





s =0

We can approximate the integrals in (75) as follows:



2π

∂ u(b, φ, λ) cos n(φ − θj )dφ Rn (tk − λ)dλ ∂λ 0 0  l +1  1 l +1 J −1  θs+1 l l k−1  tl+1   ulI+ uI ,s − ulI ,s ,s+1 − uI ,s+1 − uI ,s + uI ,s φ − θs = + · τ τ δ s =0 θ s l = 0 tl tk



· cos n(φ − θj )dφ Rn (tk − λ)dλ  tl+1 J −1 k−1  δ l +1 l Rn (tk − λ)dλ. = (u − uI ,s )[2 cos n(θs − θj ) − cos n(θs+1 − θj ) − cos n(θs−1 − θj )] τ l=0 s=0 I ,s tl We will compute the error E = max Uik,j − uki,j  ,





i,j,k

where Uik,j results.

and uki,j are the exact and numerical solutions in (ri , θj , tk ), respectively, and will be listed in the tables of numerical

Example 1. Let us consider the problem (62)–(65) with

     5  2 3 2 Γ (4 + α) 2  ( r − 1 ) exp (− r ) t ( r − 1 ) + 27 + 11r − r + tα   6 r   f (r , θ , t ) = −4r 2 (r − 1)2 16r 2 − 17r + 4 , (r , θ , t ) ∈ ΩiT ,     0, (r , θ , t ) ∈ ΩeT , and with the exact solution u(r , θ , t ) =

   (r − 1)4 exp(−r )t 3+α + r 4 , 0,

(r , θ , t ) ∈ ΩiT ,

(r , θ , t ) ∈ ΩeT ,

where a = 0.1 and b = 1. For different r step numbers I and α for T = 1, the numerical results are reported in Tables 1–4. In Figs. 1 and 2, we give the error function |U (b, θ , t ) − u(b, θ , t )| for t ∈ [0, T ], θ = 0 and different α . Example 2. Let us consider the problem (62)–(65) with

f (r , θ , t ) =

     5  2 2 3 2 Γ (4 + α)  + − 27 + 11r − r + tα ( r − 1 ) exp (− r ) t ( r − 1 )   6

+ r 2 (r − 1)2 cos θ (−15 + 66r − 63r 2 ),     0, (r , θ , t ) ∈ ΩeT ,

r

(r , θ , t ) ∈ ΩiT ,

and the exact solution u(r , θ , t ) =

 (r − 1)4 exp(−r )t 3+α + r 4 (r − 1)4 cos θ , 0, (r , θ , t ) ∈ ΩeT ,

(r , θ , t ) ∈ ΩiT ,

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

23

Table 1 Max error for α = 1, Example 1.

α

I

N =0

1

2 4 8 16

8.536162e−2 1.387306e−2 2.680616e−3 5.988132e−4

Table 2 Max error for α = 2/3, Example 1.

α

I

N =0

2/3

2 4 8 16

7.800210e−2 1.018042e−2 3.438102e−3 1.612343e−3

α

I

N =0

1/2

2 4 8 16

7.523328e−2 9.875277e−3 3.642875e−3 1.639176e−3

α

I

N =0

1/3

2 4 8 16

7.179382e−2 9.651485e−3 3.634604e−3 1.560218e−3

Table 3 Max error for α = 1/2, Example 1.

Table 4 Max error for α = 1/3, Example 1.

Table 5 Max error for α = 1, Example 2.

α

I

N =0

N =1

1

2 4 8 16

8.414797e−2 1.426294e−2 2.734573e−3 6.029104e−4

7.867256e−2 1.370433e−2 2.679262e−3 5.974965e−4

Table 6 Max error for α = 2/3, Example 2.

α

I

N =0

N =1

2/3

2 4 8 16

7.668599e−2 1.077837e−2 3.415905e−3 1.607633e−3

7.168193e−2 1.006833e−2 3.415905e−3 1.607633e−3

where a = 0.1 and b = 1. For different r step numbers I and α for T = 1, the numerical results are reported in Tables 5–8. In Figs. 3 and 4, we give the error function |U (b, θ , t ) − u(b, θ , t )| for t ∈ [0, T ], θ = π /3 and different α . 5. Conclusions In this paper, we obtained the exact and approximating artificial boundary conditions for the fractional sub-diffusion equation on an unbounded domain in two space dimensions by using Laplace transform and Fourier series expansion jointly. By some numerical examples the effectiveness of the artificial boundary conditions is tested. The original problem on unbounded domain is reduced to an initial–boundary value problem on a bounded domain with an artificial boundary by imposing one or more of the obtained artificial boundary conditions. Finally, the resulting problem is solved by the classical Crank–Nicolson method, although, any other suitable numerical method can be used similarly.

24

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

Fig. 1. The error curves with different I, number of r segments, Example 1. (left) α = 1, (right) α = 2/3.

Fig. 2. The error curves with different I, number of r segments, Example 1. (left) α = 1/2, (right) α = 1/3.

Fig. 3. The error curves with different I, number of r segments, Example 2. (left) α = 1, (right) α = 2/3. Table 7 Max error for α = 1/2, Example 2.

α

I

N =0

N =1

1/2

2 4 8 16

7.396049e−2 1.039335e−2 3.640399e−3 1.634520e−3

6.921877e−2 9.724701e−3 3.638036e−3 1.634520e−3

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

25

Fig. 4. The error curves with different I, number of r segments, Example 2. (left) α = 1/2, (right) α = 1/3. Table 8 Max error for α = 1/3, Example 2.

α

I

N =0

N =1

1/3

2 4 8 16

7.083991e−2 1.010630e−2 3.675442e−3 1.559066e−3

6.643365e−2 9.468618e−3 3.667519e−3 1.558577e−3

Acknowledgements The authors like to thank the anonymous reviewers for their comments and suggestions that have improved the paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000. R. Metzler, J. Klafter, The random walks guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (1) (2000) 1–77. R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini, P. Paradisi, Discrete random walk models for space–time fractional diffusion, Chem. Phys. 284 (1–2) (2002) 521–541. B.I. Henry, S.L. Wearne, Fractional reaction–diffusion, Physica A 276 (3–4) (2000) 448–455. M. Giona, H.E. Roman, Fractional diffusion equation for transport phenomenna in random media, Physica A 185 (1–4) (1992) 87–97. K. Diethelm, The Analysis of Fractional Differential Equations, Springer-Verlag, Berlin, 2004. C.M. Chen, F. Liu, I. Turner, V. Anh, Numerical schemes and multivariate extrapolation of a two-dimensional anomalous sub-diffusion equation, Numer. Algorithms 54 (1) (2010) 1–21. P. Zhuang, F. Liu, Finite difference approximation for two-dimensional time fractional differential diffusion equation, J. Algorithms Comput. Technol. 1 (1) (2007) 1–15. Q. Liu, Y.T. Gu, P. Zhuang, F. Liu, Y.F. Nie, An implicit RBF meshless approach for time fractional diffusion equations, Comput. Mech. 48 (1) (2011) 1–12. Y.N. Zhang, Z.Z. Sun, Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation, J. Comput. Phys. 230 (24) (2011) 8713–8728. Y.M. Wang, Maximum norm error estimates of ADI methods for a two dimensional fractional sub-diffusion equation, Adv. Math. Phys. (2013) 1–10. Article ID 293706. M. Cui, Convergence analysis of high-order compact alternating direction implicit schemes for the two dimensional time fractional diffusion equation, Numer. Algorithms 62 (3) (2013) 383–409. X. Yang, H. Zhang, D. Xu, Orthogonal spline collocation method for the two-dimensional fractional sub-diffusion equation, J. Comput. Phys. 256 (2014) 824–837. L. Yan, F. Yang, A kansa-type MFS scheme for two dimensional time fractional diffusion equations, Eng. Anal. Bound. Elem. 37 (2013) 1426–1435. H. Jiang, F. Liu, I. Turner, 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 (2012) 1117–1127. H. Jiang, F. Liu, I. Turner, K. Burrage, Analytical solutions for the generalized multi-term time-fractional diffusion-wave/diffusion equation in a finite domain, Comput. Math. Appl. 64 (2012) 3377–3388. H. Jiang, F. Liu, M.M. Meerschaert, R. McGough, Fundamental solutions for the multi-term modified power law wave equations in a finite domain, Electron. J. Math. Anal. Appl. 1 (7) (2013) 1–12. F. Huang, F. Liu, The time fractional diffusion equation and the advection–dispersion equation, ANZIAM J. 46 (2005) 317–330. F. Huang, F. Liu, The fundamental solution of the space–time fractional advection–dispersion equation, Appl. Math. Comput. 18 (1–2) (2005) 339–350. F. Huang, F. Liu, The space–time fractional diffusion equation with Caputo derivatives, Appl. Math. Comput. 19 (1–2) (2005) 179–190. D. Qikui, T. Minxia, Exact and approximate artificial boundary conditions for the hyperbolic problems in unbounded domain, Appl. Math. Comput. 169 (1) (2005) 544–562. B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (139) (1977) 629–651. T. Hagstrom, S.I. Hariharan, D. Thompson, High-order radiation boundary conditions for the convective wave equation in exterior domains, SIAM J. Sci. Comput. 25 (3) (2003) 1088–1101. X. Wu, Z.Z. Sun, Convergence of difference scheme for heat equation in unbounded domains using artificial boundary conditions, Appl. Numer. Math. 50 (2) (2004) 261–277.

26

R. Ghaffari, S.M. Hosseini / Computers and Mathematics with Applications 68 (2014) 13–26

[26] H. Han, Z. Huang, Exact and approximating boundary conditions for the parabolic problems on unbounded domains, Comput. Math. Appl. 44 (5–6) (2002) 655–666. [27] H.D. Han, Z.Y. Huang, A class of artificial boundary conditions for heat equation in unbounded domains, Comput. Math. Appl. 43 (6–7) (2002) 889–900. [28] H. Han, X. Wu, Artificial Boundary Method, Tsinghua Univ. Press, 2013. [29] G.H. Gao, Z.Z. Sun, Y.N. Zhang, A finite difference scheme for fractional sub-diffusion equations on an unbounded domain using artificial boundary conditions, J. Comput. Phys. 231 (7) (2012) 2865–2879. [30] G.H. Gao, Z.Z. Sun, The finite difference approximation for a class of fractional sub-diffusion equations on a space unbounded domain, J. Comput. Phys. 236 (2013) 443–460. [31] S.M. Hosseini, R. Ghaffari, Polynomial and nonpolynomial spline methods for fractional sub-diffusion equations, Appl. Math. Model. (2014) http://dx.doi.org/10.1016/j.apm.2013.11.062. [32] L. Debnath, Integral Transforms and their Applications, CRC Press, 1995. [33] M. Abramovitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, 1964. [34] H.S. Carslaw, J.C. Jaeqer, Operational Methods in Applied Mathematics, Oxford Univ. Press, New York, 1963. [35] G.N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge Univ. Press, Cambridge, 1966. [36] B.G. Kornev, Bessel Functions and their Applications, CRC Press, New York, 1798. [37] A. Gray, G.B. Mathews, A Treatise on Bessel Functions and their Applications to Physics, London Macmillan, New York, 1895. [38] A.A.M. Cuyt, V. Petersen, B. Verdonk, H. Waadeland, W.B. Jones, Handbook on Continued Fractions for Special Functions, Springer, 2008. [39] F.B. Hildebrand, Advanced Calculus for Applications, Prentice Hall, 1948. [40] G. Doetsch, Introduction to the Theory and Application of the Laplace Transformation, Springer-Verlag, Berlin, Heidelberg, New York, 1974. [41] M.K. Kerimove, S.L. Skorokhodov, Calculation of the complex zeros of the modified bessel functions of the second kind and its derivatives, Comput. Math. Math. Phys. 24 (4) (1984) 115–123. [42] R. Parnes, Complex zeros of the modified bessel functions Kn (z ), Math. Comp. 26 (120) (1972) 949–953. [43] D.G. Randall, Supersonic flow past quasicylindrical bodies of almost circular crosssection, Reports and Memoranda—Aeronautical Research Council, 3067, 1958. [44] H.F. Weinberger, A First Course in Partial Differential Equations with Complex Variables and Transform Methods, University of Minnesota, New York, 1965. [45] H.S. Carslaw, J.C. Jaeqer, Conduction of Heat in Solids, Clarendon Press, Oxford, 1959.