Numerical solution for a sub-diffusion equation with a smooth kernel

QUT Digital Repository: http://eprints.qut.edu.au/ Liu, F. and Yang, C. and Burrage, K. (2009) Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term. Journal of Computational and Applied Mathematics, 231(1). pp. 160176. © Copyr igh t 2 0 0 9 Else vie r Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term ⋆ F. Liu a,c,∗ , C. Yang b , K. Burrage c,d a School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld. 4001, Australia b School of Mathematical Sciences, Xiamen University, Xiamen 361005, China c IMB, University of Queensland, Qld. 4072, Australia. d COMLAB and OCISB, Oxford University, OX1 3LB, UK. Abstract In this paper, we consider a modified anomalous subdiffusion equation with a nonlinear source term for describing processes that become less anomalous as time progresses by the inclusion of a second fractional time derivative acting on the diffusion term. A new implicit difference method is constructed. The stability and convergence are discussed using a new energy method. Finally, some numerical examples are given. The numerical results demonstrate the effectiveness of theoretical analysis. Key words: Implicit difference method, modified anomalous subdiffusion equation, nonlinear source terms, energy method, stability and convergence. MSC(2000) 26A33, 34K28, 65M12, 60J70 1 Introduction In recent years, it has been reported that, in numerous physical and biological systems many diffusion rates of species cannot be characterized by the single ⋆ This research has been supported by the National Natural Science Foundation of China grant 10271098. ∗ Corresponding author. Email address: f.liu@qut.edu.au or fwliu@xmu.edu.cn (F. Liu). Preprint submitted to Elsevier Science 9 November 2008 parameter of the diffusion constant [40]. Instead, the (anomalous) diffusion is characterized by a scaling parameter γ as well as a diffusion constant K, and the mean square displacement of diffusing species hx2 (t)i scales as a nonlinear power-law in time, i.e., 2Kγ γ t , t → ∞, Γ(1 + γ) where γ (with 0 < γ < 1) is the anomalous diffusion exponent and Kγ is the generalized diffusion coefficient. Ordinary (or Brownian) diffusion corresponds to γ = 1 with K1 = D (the ordinary diffusion coefficient). For example, single particle tracking experiments and photo-bleaching recovery experiments have revealed sub-diffusion (0 < γ < 1) of proteins and lipids in a variety of cell membranes [2,7–9,35,37,38]. Anomalous subdiffusion has also been observed in neural cell adhesion molecules [36]. Indeed anomalous subdiffusion (the case with 0 < γ < 1) is generic in media with obstacles [24,25] or binding sites [26]. For anomalous subdiffusive random walks, the continuum description via the ordinary diffusion equation is replaced by the fractional diffusion equation. It has been suggested that the probability density function (pdf) u(x, t) that describes anomalous subdiffusion particles follows the anomalous sub-diffusion equation [21]: hx2 (t)i ∼ ∂u ∂2u ∂ 1−γ = 1−γ Kγ 2 , 0 ≤ x ≤ a, 0 < t ≤ T, ∂t ∂t ∂x " # (1) where u(x, t) is the probability density that the particle that started at 0 at time 0 is at x at time t, ∂ 1−γ u/∂t1−γ denotes the Riemann-Liouville fractional derivative of order 1 − γ defined by 1 ∂ Z t u(x, η) ∂ 1−γ 1−γ u(x, t) = u(x, t) = D dη, 0 t ∂t1−γ Γ(γ) ∂t 0 (t − η)1−γ (2) with 0 ≤ γ ≤ 1. For γ = 1 one recovers the identity operator, and for γ = 0 the ordinary first-order derivative. Yuste and Lindenberg considered a combination of these two phenomena and proposed to solve the A+A reaction-subdiffusion problem in one dimension [43]. Further, Yuste, Acedo and Lidenberg [44] proposed the A + B reactionsubdiffusion equations: ∂ ∂2 a(x, t) = Kγ 0 Dt1−γ 2 a(x, t) − Rγ (x, t), ∂t ∂x ∂ ∂2 b(x, t) = Kγ 0 Dt1−γ 2 b(x, t) − Rγ (x, t). ∂t ∂x The reaction term has many different forms. Seki, Wojcik and Tachiya [34] 2 proposed a reaction-subdiffusion equation which at long times corresponds to choosing a reaction term of the form Rγ (x, t) = k0 Dt1−γ a(x, t)b(x, t). Tan et al. [41] and Chen et al. [5] considered Stokes’ first problem for a heated generalized second grade fluid with fractional derivative with a nonhomogeneous forcing term: ∂ 2 u(x, t) ∂ 2 u(x, t) ∂u(x, t) + κ =0 Dt1−γ κ1 + f (x, t). 2 ∂t ∂x2 ∂x2 " # (3) There are numerous approaches to modeling anomalous diffusive behaviour, such as, continuous time random walks, Monte Carlo simulations [25], Langevin equations and fractional diffusion equations [21]. The fractional diffusion equation is characterised by the presence of either a fractional temporal derivative or fractional spatial derivative or both (time-fractional diffusion equations were introduced by Zaslavsky [46], and also references [21,22] for a recent review). Other fractional variants are the fractional Fokker-Planck equation for anomalous diffusion due to an external force and fractional reaction-diffusion equations [6,10,32,33,44] for reactions where the products and reactants diffuse anomalously. These equations involve only a single temporal fractional derivative acting on the diffusion term. Recently a model has been proposed [3,12,39,40] for describing processes that become less anomalous as time progresses by the inclusion of a secondary fractional time derivative acting on a diffusion operator, Lx = K∂ 2 /∂x2 , ∂u(x, t) ∂ 1−α ∂ 1−β = A 1−α + B 1−β ∂t ∂t ∂t ! Lx u(x, t), (4) where 0 < α < β 6 1, K is diffusion coefficient, A and B are positive constants. The subdiffusive motion is characterized by an asymptotic longtime behavior of the mean square displacement of the form hx2 (t)i = 2A 2B tα + tβ . Γ(1 + α) Γ(1 + β) (5) A possible application of this equation is in econophysics where there is an increasing interest in modelling using continuous time random walks [16,23,27– 31]. In particular the crossover between more and less anomalous behaviour has been observed in the volatility of some share prices [17–19]. In this paper, we consider the following modified anomalous subdiffusion equation with a nonlinear source term: ∂ 1−α ∂ 1−β ∂u(x, t) = A 1−α + B 1−β ∂t ∂t ∂t !" 3 ∂ 2 u(x, t) + f (x, t) + g(u, x, t). ∂x2 # (6) Recently, many researchers have proposed various numerical methods to solve the space or time fractional partial differential equations. Liu, Anh and Turner [13] proposed a computationally effective method of lines for the space fractional partial differential equation. They transformed the space fractional partial differential equation into a system of ordinary differential equations that was then solved using backward differentiation formulas. Meerschaert and Tadjeran [20] developed finite difference approximations for fractional advectiondispersion flow equations. Roop [14] investigated the computational aspects of the Galerkin approximation using continuous piecewise polynomial basis functions on a regular triangulation of a bounded domain in R2 . Liu et al. [15] also investigated the stability and convergence of difference methods for the space-time fractional advection-diffusion equation. Yu et al. [42] developed a reliable algorithm of the Adomian decomposition method to solve the linear and nonlinear space-time fractional reaction-diffusion equations in the form of a rapidly convergent series with easily computable components. Yuste and Acedo [45] proposed an explicit finite difference method and a new Von Neumann-type stability analysis for the anomalous sub-diffusion equation (1). However, they did not give a convergence analysis and pointed out the difficulty of this task when implicit methods are considered. Langlands and Henry [11] also investigated this problem and proposed an implicit numerical scheme (L1 approximation), and discussed the accuracy and stability of this scheme. However, the global accuracy of the implicit numerical scheme has not been derived and it seems that the unconditional stability for all γ in the range 0 < γ ≤ 1 has not been established. Recently, Chen and Liu et al. [4] presented a Fourier method for the anomalous sub-diffusion equation (1), and they gave the stability analysis and the global accuracy analysis of the difference approximation scheme. Zhuang and Liu et al. [47] also proposed new solution and analytical techniques of implicit numerical methods for the anomalous subdiffusion equation (1). Chen and Liu et al. [5] proposed implicit and explicit numerical approximation schemes for the Stokes’ first problem for a heated generalized second grade fluid with fractional derivative (3). The stability and convergence of the numerical schemes are discussed using a Fourier method. A Richardson extrapolation technique for improving the order of convergence of the implicit scheme is presented. However, effective numerical methods and error analysis for the modified anomalous subdiffusion equation with a nonlinear source term are still in their infancy and are open problems. The main purpose of this paper is to solve and analyze this problem by introducing an implicit difference method and new analytical techniques. The structure of the remainder of this paper as follows. In Section 2, an implicit numerical method for the modified anomalous subdiffusion equation with a nonlinear source term is proposed. In Sections 3 and 4, the stability and convergence of the implicit numerical method are discussed, respectively. Finally, some numerical results for the modified anomalous subdiffusion equation with 4 a nonlinear source term are given. 2 An implicit numerical method for the modified anomalous subdiffusion equation In this paper, we consider the function describing processes that get less anomalous in the course of time, namely ∂u(x, t) ∂ 2 u(x, t) ∂ 1−α ∂ 1−β = A 1−α + B 1−β + f (x, t) + g(u, x, t), ∂t ∂t ∂ ∂x2 0 ≤ x ≤ L, 0 ≤ t ≤ T, !" # (7) with initial and boundary conditions: u(x, 0) = φ(x), 0 ≤ x ≤ L, u(0, t) = ϕ1 (t), u(L, t) = ϕ2 (t), 0 ≤ t ≤ T, (8) (9) where 0 ≤ α ≤ β ≤ 1, u(x, t) is the probability density of the particle that started at time 0 is at x at time t. In Eq. (7), the operators ∂ 1−α /∂t1−α and ∂ 1−β /∂t1−β denote the Riemann-Liouville fractional derivative operators. Let Ω = [0, L] × [0, T ]. We define the function space G(Ω) = {w(x, t)| ∂5w ∂ 2w 2 ∈ C (Ω) and ∈ C(Ω)}. ∂x2 ∂x4 ∂t In this paper, we suppose the continuous problem (7-9) has a smooth solution u(x, t) ∈ G(Ω). Now we construct an implicit difference method using a new solution technique. We define tk = kτ, k = 0, 1, . . . , n; xi = ih, i = 0, 1, . . . , m, where τ = T /n and h = L/m are space and time step sizes, respectively. We introduce the following notations: ∂ 2 u(x, t) , ∂x2 δx2 u(x, t) = u(x + h, t) − 2u(x, t) + u(x − h, t), δ 2 u(x, t) Lh v(x, t) = x 2 . h Lu(x, t) = Integrating both sides of Eq. (7) from tk to tk+1 , we have 5 (10) (11) (12) u(xi , tk+1 ) − u(xi , tk ) A Z tk+1 Lu(xi , η) + f (xi , η) A Z tk Lu(xi , η) + f (xi , η) dη − dη = Γ(α) 0 (tk+1 − η)1−α Γ(α) 0 (tk − η)1−α B Z tk+1 Lu(xi , η) + f (xi , η) B Z tk Lu(xi , η) + f (xi , η) + dη − dη Γ(β) 0 (tk+1 − η)1−β Γ(β) 0 (tk − η)1−β + Z tk+1 tk g(u(xi , η), xi , η)dη. Hence u(xi , tk+1 ) = u(xi , tk ) A Z τ Lu(xi , η) + f (xi , η) dη + + Γ(α) 0 (tk+1 − η)1−α B Z τ Lu(xi , η) + f (xi , η) + dη + Γ(β) 0 (tk+1 − η)1−β + Z tk+1 tk A Z tk Lv(xi , η) + p(xi , η) dη Γ(α) 0 (tk − η)1−α B Z tk Lv(xi , η) + p(xi , η) dη Γ(β) 0 (tk − η)1−β g(u(xi , η), xi , η)dη, (13) where v(x, t) = u(x, t + τ ) − u(x, t), p(x, t) = f (x, t + τ ) − f (x, t). Let B Z τ Lu(xi , η) + f (xi , η) A Z τ Lu(xi , η) + f (xi , η) dη + dη, Γ(α) 0 (tk+1 − η)1−α Γ(β) 0 (tk+1 − η)1−β A Z tk Lv(xi , η) + p(xi , η) B Z tk Lv(xi , η) + p(xi , η) I2 = dη + dη, Γ(α) 0 (tk − η)1−α Γ(β) 0 (tk − η)1−β I1 = I3 = Z tk+1 tk g(u(xi , η), xi , η)dη. Eq. (13) can be written as u(xi , tk+1 ) = u(xi , tk ) + I1 + I2 + I3 . For I1 , we can get approximation as below: 6 A Z τ Lu(xi , τ ) + f (xi , τ ) I1 = dη Γ(α) 0 (tk+1 − η)1−α B Z τ Lu(xi , τ ) + f (xi , τ ) + dη + R11 Γ(β) 0 (tk+1 − η)1−β Aτ α [(k + 1)α − k α ] [Lu(xi , τ ) + f (xi , τ )] = Γ(α + 1) i Bτ β h + (k + 1)β − k β [Lu(xi , τ ) + f (xi , τ )] + R11 Γ(β + 1)   1 2 Aτ α ak 2 δx u(xi , τ ) + f (xi , τ ) + R12 = Γ(α + 1) h   β Bτ 1 2 + bk 2 δx u(xi , τ ) + f (xi , τ ) + R12 + R11 Γ(β + 1) h   α Aτ 1 2 ak 2 δx u(xi , τ ) + f (xi , τ ) = Γ(α + 1) h   β Bτ 1 + bk 2 δx2 u(xi , τ ) + f (xi , τ ) + R1 Γ(β + 1) h # "  α Bτ β bk 1 2 Aτ ak + δ u(xi , τ ) + f (xi , τ ) + R1 , = Γ(α + 1) Γ(β + 1) h2 x where ak = (k + 1)α − k α , bk = (k + 1)β − k β , R1 = R11 + Aτ α ak Bτ β bk R12 + R12 . Γ(α + 1) Γ(β + 1) A Z τ Lu(xi , η) − Lu(xi , τ ) + f (xi , η) − f (xi , τ ) dη Γ(α) 0 (tk+1 − η)1−α B Z τ Lu(xi , η) − Lu(xi , τ ) + f (xi , η) − f (xi , τ ) + dη, Γ(β) 0 (tk+1 − η)1−β 1 ∂ 2 u(xi , τ ) − 2 δx2 u(xi , τ ). R12 = 2 ∂x h R11 = Given that Lu(xi , η) + f (xi , η) # " ∂ 2 u(xi , τ ) ∂ 3 u(xi , ξ1 ) ∂f (xi , ξ2 ) (η − τ ), = + f (xi , τ ) + + ∂x2 ∂x2 ∂t ∂t where 0 ≤ τ ≤ ξ1 ≤ η; 0 ≤ τ ≤ ξ2 ≤ η, we obtain 7 (14) Z τ Z τ B 1 1 A C1 τ dη + C2 τ dη |R11 | ≤ 1−α Γ(α) Γ(β) 0 (tk+1 − η) 0 (tk+1 − η)1−β AC1 τ 1+α BC2 τ 1+β ≤ ak + bk . (15) Γ(α + 1) Γ(β + 1) Again, it is apparent that |R12 | ≤ C3 h2 . So, we have |R1 | ≤ BC2 τ 1+α Aτ α Bτ β AC1 τ 1+α ak + bk + ak C3 h2 + bk C3 h2 Γ(α + 1) Γ(β + 1) Γ(α + 1) Γ(β + 1)    ≤ Cak τ α τ + h2 + Cbk τ β τ + h2  ≤ C ak τ α + b k τ β   τ + h2 .  (16) For I2 , we get the following approximation I2 = X Z tj+1 Lv(xi , η) + p(xi , η) A k−1 dη Γ(α) j=0 tj (tk − η)1−α X Z tj+1 Lv(xi , η) + p(xi , η) B k−1 + dη Γ(β) j=0 tj (tk − η)1−β = X Z tj+1 Lv(xi , tj+1 ) + p (xi , tj+1 ) A k−1 dη Γ(α) j=0 tj (tk − η)1−α X Z tj+1 Lv(xi , tj+1 ) + p (xi , tj+1 ) B k−1 + dη + R21 Γ(β) j=0 tj (tk − η)1−β = X Aτ α k−1 ak−j−1 [Lh v(xi , tj+1 ) + p (xi , tj+1 ) + R22 ] Γ(α + 1) j=0 +  = X Bτ β k−1 bk−j−1 [Lh v(xi , tj+1 ) + p (xi , tj+1 ) + R22 ] + R21 Γ(β + 1) j=0 α k−1 X β k−1 X  Bτ Aτ ak−j−1 + bk−j−1  Γ(α + 1) j=0 Γ(β + 1) j=0 × [Lh v(xi , tj+1 ) + p (xi , tj+1 )] + R2 , where 8 (17) R2 = R21 + X X Bτ β k−1 Aτ α k−1 ak−j−1 R22 + bk−j−1 R22 , Γ(α + 1) j=0 Γ(β + 1) j=0 X Z tj+1 Lv(xi , η) − Lv(xi , tj+1 ) + p(xi , η) − p(xi , tj+1 ) A k−1 dη R21 = Γ(α) j=0 tj (tk+1 − η)1−α + X Z tj+1 Lv(xi , η) − Lv(xi , tj+1 ) + p(xi , η) − p(xi , tj+1 ) B k−1 dη, Γ(β) j=0 tj (tk+1 − η)1−β R22 = Lv(xi , tj+1 ) − Lh v(xi , tj+1 ). (18) Because Lv(xi , η) + p(xi , η) # " ∂ 3 v(xi , η1 ) ∂p(xi , η2 ) ∂ 2 v(xi , tj+1 ) (η − tj+1 ) + p (xi , tj+1 ) + + = ∂x2 ∂x2 ∂t ∂t # " ∂ 2 v(xi , tj+1 ) ∂ 4 v(xi , η˜1 ) ∂ 2 p(xi , η˜2 ) = τ (η − tj+1 ), + p (xi , tj+1 ) + + ∂x2 ∂x2 ∂t2 ∂t2 where η ≤ η1 ≤ tj+1 , we have η1 ≤ η˜1 ≤ η1 + τ, η ≤ η2 ≤ tj+1 , η2 ≤ η˜2 ≤ η2 + τ , X Z tj+1 X Z tj+1 1 1 B 2 k−1 A 2 k−1 τ dη + τ dη |R21 | ≤ 1−α Γ(α) j=0 tj (tk − η) Γ(β) j=0 tj (tk − η)1−β ! A B ≤ kα + kβ τ 2. Γ(α + 1) Γ(β + 1) And, using Taylor’s formula, we can obtain h2 ∂ 4 v(ξ2 , tj+1 ) 12 " ∂x4 # h2 ∂ 4 u(ξ2 , tj+2 ) ∂ 4 u(ξ2 , tj+1 ) = Lh v(xi , tj+1 ) + − 12 ∂x4 ∂x4 h2 τ ∂ 5 u(ξ2 , η̃3 ) . = Lh v(xi , tj+1 ) + 12 ∂x4 ∂t Lv(xi , tj+1 ) = Lh v(xi , tj+1 ) + Hence, we have |R22 | ≤ |Lv(xi , tj+1 ) − Lh v(xi , tj+1 )| ≤ Cτ h2 . 9 (19) X X Aτ α k−1 Bτ β k−1 |R2 | ≤ R21 + ak−j−1 R22 + bk−j−1 R22 Γ(α + 1) j=0 Γ(β + 1) j=0   Bτ β+1 β  Aτ α+1 α  k τ + h2 + k τ + h2 Γ(α + 1) Γ(β + 1)     Aτ Bτ T α τ + h2 + T β τ + h2 ≤ Γ(α + 1) Γ(β + 1) ≤   ≤ Cτ τ + h2 . (20) For I3 , we can get the approximation as below: I3 = = Z tk+1 tk g(u(xi , η), xi , η)dη τ [g(u(xi , tk+1 ), xi , tk+1 ) + g(u(xi , tk ), xi , tk )] + O(τ 2 ). 2 (21) From the above result, we obtain u(xi , tk+1 ) Bτ β bk Aτ α ak = u(xi , tk ) + + Γ(α + 1) Γ(β + 1) " # 1 2 δ u(xi , τ ) + f (xi , τ ) h2 x    X X Bτ β k−1 1 Aτ α k−1 ak−j−1 + bk−j−1  2 δx2 v(xi , tj+1 ) + p (xi , tj+1 ) + Γ(α + 1) j=0 Γ(β + 1) j=0 h τ + [g(u(xi , tk+1 ), xi , tk+1 ) + g(u(xi , tk ), xi , tk )] + Rik+1 2 # "  1 2 Bτ β bk Aτ α ak + δ u(xi , τ ) + f (xi , τ ) = u(xi , tk ) + Γ(α + 1) Γ(β + 1) h2 x + k−1 X j=0 + " Bτ β Aτ α ak−j−1 + bk−j−1 Γ(α + 1) Γ(β + 1) # 1 2 δ v(xi , tj+1 ) + p (xi , tj+1 ) h2 x τ [g(u(xi , tk+1 ), xi , tk+1 ) + g(u(xi , tk ), xi , tk )] + Rik+1 . 2 Let Aτ α Bτ β = r1 , = r2 . Γ(α + 1) Γ(β + 1) We have 10    (22) u(xi , tk+1 ) h i = u(xi , tk ) + (r1 ak + r2 bk ) δx2 u(xi , τ )/h2 + f (xi , τ ) + k−1 X (r1 ak−j−1 + r2 bk−j−1 ) × j=0 h  δx2 u(xi , tj+2 ) − δx2 u(xi , tj+1 ) /h2 +f (xi , tj+2 ) − f (xi , tj+1 )] τ + [g (u(xi , tk+1 ), xi , tk+1 ) + g(u(xi , tk ), xi , tk )] + Rik+1 , 2 (23) where  Rik+1 ≤ C ak τ α + bk τ β     τ + h2 + Cτ τ + h2 ≤ C̃ ak τ α + bk τ β + τ   τ + h2 .  (24) From [47] we have following lemma. Lemma 1: In (24), the coefficient ak , bk (k = 0, 1, 2 . . .) satisfy: (1) a0 = 1, ak > 0, b0 = 1, bk > 0, k = 0, 1, 2, . . . ; (2) ak > ak+1 , bk > bk+1 , k = 0, 1, 2, . . . ; (3) there exists a positive constant C > 0, such that τ ≤ Cak τ α , τ ≤ Cbk τ β , k = 1, 2, . . . . Let u = (u1 , u2 , . . . , um−1 )T , v = (v1 , v2 , . . . , vm−1 )T . We define (u, v) = m−1 X ui vi h, ku2 k = i=1 q (u, u) = m−1 X i=1 ! 12 u2i h . We have v u m−1 u X k k k k T k |Rik |2 . R = (R1 , R2 , . . . , Rm−1 ) , kR k2 = th i=1 We obtain Rk uki 2  ≤ C ak τ α + b k τ β   τ + h2 . We denote for the numerical approximation to u(xi , tk ), fik = f (xi , tk ), k gi = (u(xi , tk ), xi , tk ) and write δx2 uki = uki+1 − 2uki + uki−1 . We obtain the following implicit difference scheme: 11 h uk+1 = uki + (r1 ak + r2 bk ) δx2 u1i /h2 + fi1 i + k−1 X h i (r1 ak−j−1 + r2 bk−j−1 ) (δx2 uj+2 − δx2 uj+1 )/h2 + fij+2 − fij+1 i i j=0  τ + gik+1 + gik . 2 i (25) The implicit numerical method can be rewritten as the following form: uk+1 = uki + r1 a0 (δx2 uk+1 /h2 + fik+1 ) + r1 i i k−1 X (aj+1 − aj )(δx2 uik−j /h2 + fik−j ) j=0 +r2 b0 (δx2 uk+1 /h2 + fik+1 ) + r2 i k−1 X (bj+1 − bj )(δx2 uik−j /h2 + fik−j ) j=0 τ + (gik+1 + gik ) 2 = uki + (r1 + r2 )(δx2 uk+1 /h2 + fik+1 ) i + k−1 X [r1 (aj+1 − aj ) + r2 (bj+1 − bj )](δx2 uik−j /h2 + fik−j ) j=0 τ + (gik+1 + gik ). 2 (26) The initial and boundary conditions are u0i = φ(ih), i = 0, 1, 2, . . . , m; uk0 = ϕ1 (kτ ), ukm = ϕ2 (kτ ), k = 1, 2, . . . , n. (27) (28) We simplify the equation as  uk+1 − uki = (r1 + r2 ) δx2 uk+1 /h2 + fik+1 i i + k−1 X   [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] δx2 uik−j /h2 + fik−j j=0  τ + gik+1 + gik . 2  (29) We sum from u1i − u0i to uk+1 − uki , and use the definition of aj , bj to obtain i 12 uk+1 − u0i = i k−1 X  (r1 aj+1 + r2 bj+1 ) δx2 uik−j /h2 + fik−j j=−1 + =   X j τ k+1 gi + gij−1 2 j=1 k X   (r1 aj + r2 bj ) δx2 uk+1−j /h2 + fik+1−j + τ i j=0  τ + gi0 + gik+1 . 2 k X gij j=1 (30) So that, we have: uki = u0i + k−1 X (r1 aj + r2 bj ) j=0  δx2 uik−j /h2 + fik−j  +τ k−1 X τ gij + (gi0 + gik ), (31) 2 j=1 where Aτ α Bτ β , r2 = , Γ(α + 1) Γ(β + 1) ak = (k + 1)α − k α , bk = (k + 1)β − k β , i = 1, 2, . . . , m − 1, k = 1, 2, . . . , n. The initial and boundary conditions are (27),(28). δx2 uki = uki+1 − 2uki + uki−1 , r1 = 3 Stability of the implicit numerical method We give a stability analysis as follows. We suppose that ũki , i = 0, 1, . . . , m; k = 0, 1, . . . , n is the approximate solution of Eq.(7), g̃ik denotes g (ũ(xi , tk ), xi , tk ), the error εki = uki − ũki , satisfies εk0 = 0, εkm = 0, k = 1, 2, . . . , n and  εk+1 = εki + (r1 + r2 ) × δx2 εk+1 /h2 i i + k−1 X   [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × δx2 εik−j /h2 j=0  τ + gik+1 − g̃ik+1 + gik − g̃ik . 2  (32) We also suppose that the function g(u(x, t), x, t) satisfies the Lipschitz condition, |g(u1 , x, t) − g(u2 , x, t)| ≤ L|u1 − u2 |, ∀u1 , u2 , so that gik − g̃ik ≤ L εki . 13 We can easily prove the following result: Lemma 2: Let ∆vi = vi+1 −vi , ∆wi = wi+1 −wi , δ 2 vi = vi+1 −2vi +vi−1 , δ 2 wi = wi+1 − 2wi + wi−1 . If v0 = wm = 0, then   δ 2 v, w = −v1 w1 h − (∆v, ∆w) where  δ 2 v = δ 2 v1 , δ 2 v2 , · · · , δ 2 vm−1 T , ∆v = (∆v1 , ∆v2 , · · · , ∆vm−1 )T , ∆w = (∆w1 , ∆w2 , · · · , ∆wm−1 )T . Proof:   δ 2 v, w = m−1 X δ 2 v i wi h i=1 m−1 X =h =h i=1 m−1 X (vi+1 − vi ) wi − h (vi+1 − vi ) wi − h m−1 X i=1 m−1 X (vi − vi−1 ) wi (vi+1 − vi ) wi+1 i=1 i=1 +h(vm − vm−1 )wm − h(v1 − v0 )w1 = −v1 w1 h − (∆v, ∆w) . Now let E k = (εk1 , εk2 , . . . , εkm−1 )T . Multiplying Eq.(32) by hεk+1 , summing for i i from 1 to m − 1 , we obtain: E k+1 2 2   = E k+1 , E k + (r1 + r2 ) × + k−1 X h  δx2 E k+1 , E k+1 /h2 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × j=0  X τ m−1 + gik+1 − g̃ik+1 + gik − g̃ik εk+1 h i 2 i=1  = E − k+1 k−1 X ,E k  − (r1 + r2 )  2 εk+1 1 h+ h i  δx2 E k−j , E k+1 /h2 2 ∆x E k+1 2  2 /h i  [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] j=0 × h  k−j , ∆x E k+1 ε1k−j εk+1 1 h + ∆x E  /h2  X τ m−1 gik+1 − g̃ik+1 + gik − g̃ik εk+1 h. + i 2 i=1 14 i (33) According to the Schwarz inequality, we have:  2  2| ∆x E k−j , ∆x E k+1 | ≤ ∆x E k−j  k−j 2|ε1k−j εk+1 1 | ≤ ε1 2 2  + ∆x E k+1 + εk+1 1 2 2 , 2 . Because k−1 X [(aj+1 − aj )] = ak − a0 = ak − 1 and ak > 0, j=0 we have E k+1 1 ≤ 2 − ×  2 2 k+1 2 E 2 + 2 Ek 2  − (r1 + r2 ) ×  X 1 k−1 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] 2 j=0  ε1k−j 2  h + εk+1 1 2 h + ∆E k−j 2 2 2 εk+1 i + 2 ∆E k+1 2 + ∆E K+1 2 2  /h2  2 /h   2 2 2 τL τL + E k+1 + E k+1 + E k 2 2 2 2 4   2 2 1 ≤ E k+1 + E k 2 2 2    2 (1 + ak )r1 + (1 + bk )r2 2 k+1 2 × /h h + εk+1 ∆E − 1 2 2    k−1   2 2 1X /h2 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × − ε1k−j + ∆E k−j 2 2 j=0 + ≤ 1 2 3τ L E k+1 4  E k+1   2 τL Ek 2 2 4    2 2 r + r2  k+1 2 1 k+1 2 k 2 ε1 h + ∆E + E /h − 2 2 2 2 2 + X 1 k−1 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × − 2 j=0 + 3τ L E k+1 4 2 2 + τL Ek 4 2 2 . Thus 15  2 ε1k−j h+ 2 ∆E k−j 2  2 /h  (34) 2 E k+1 2 Ek 2 ≤ + 2 + + k X  [r1 aj + r2 bj ] j=0 k−1 X [r1 aj + r2 bj ] j=0 3τ L E k+1 2 2 2 +  τL Ek 2 2 2 εk+1−j 1 2 ε1k−j 2 2 h + ∆E k+1−j h+ 2 ∆E k−j 2  2 2 /h  /h2   . (35) We define the energy norm 2 Ek E = 2 Ek 2 Suppose that τ < 2 , 3L  + k−1 X (r1 aj + r2 bj ) j=0  2 ε1k−j 2 ∆E k−j 2 h+  2 /h  . then we have 1− 3τ L 2  E k+1  2 ≤ 1+ E τL 2  Ek 2 . E Thus, we have 2 Ek E ≤ 1 + 12 τ L 1 − 32 τ L !k−1 2 E1 E 1 + 21 τ L 1 − 23 τ L ≤ !n E1 2 E . Note that n = T /τ , and lim n→∞ 1 + 12 Lτ 1 − 23 Lτ !n 1+ 1− = n→∞ lim 1 LT 2n 3 LT 2n !n 1 e 2 LT = − 23 LT e = e2LT . Hence, there is a positive constant C1 > 0, such that 1 + 21 Lτ 1 − 32 Lτ thereby Ek 2 2 ≤ Ek !n 2 E ≤ C1 , ≤ C1 E 1 2 E . We have   ε1i = ε0i + (r1 + r2 ) × δx2 ε1i /h2 + 2 E1 2 ≤ 2 E0 2 − (r1 + r2 ) × 3 + τ L E1 2 2 2  2 ε11 1 + τ L E0 2 2 2 . 16 h+  τ 1 gi − g̃i1 + gi0 − g̃i0 2 2 ∆E 1 2  2 /h  (36) Therefore, 2 E1 E 1 + 21 τ L 1 − 32 τ L ≤ ! E0 2 2 . Hence, we have E1 2 2 ≤ C E0 2 2 , thus giving a stability result, i.e., kE1 k2E = kE1 k22 + r1 k X bj {|ε11−j |2 h + k∆x E1−j k22 } ≤ kE0 k22 , (37) j=0 so that kEk+1 k22 ≤ kE0 k22 . Hence, the following theorem of stability is obtained. Theorem 1. Assuming τ < defined by (25) is stable. 4 2 , 3L the fractional implicit numerical method Convergence of the implicit numerical method We let u(xi , tk ), i = 0, 1, . . . , m, k = 0, 1, . . . , n be the exact solution of equation (7) at mesh point (xi , tk ). Define ηik = u(xi , tk ) − uki , i = 0, 1, . . . , m, k =  k 0, 1, . . . , n, Y k = η1k , . . . , ηm−1  T . We get ηik+1 = ηik + (r1 + r2 ) × δx2 ηik+1 /h2 + k−1 X   [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × δx2 ηik−j /h2 j=0  i τh + g (u (xi , tk+1 ) , xi , tk+1 ) − gik+1 + g (u (xi , tk ) , xi , tk ) − gik 2 +Rik+1 , (38) k where ηi0 = η0k = ηm = 0, i = 0, 1, . . . , m. Multiplying by hη k+1 , and summing for i from 1 to m − 1, we obtain: 17 2 Y k+1  = Y + 2 k+1 k−1 X  , Y k + (r1 + r2 ) × h  δx2 Y k+1 , Y k+1 /h2 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × j=0 +   δx2 Y k−j , Y k+1 /h2    X τ m−1 k+1 k+1 gik+1 − g̃ik+1 + gik − g̃ik εk+1 h + R , Y i 2 i=1  = Y k+1 , Y k − (r1 + r2 ) × − h i k−1 X  η1k+1 2 2 h + ∆x Y k+1 2  /h2 i  [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] j=0 × + h  η1k−j η1k+1 h + ∆x Y k−j , ∆x Y k+1  /h2 i    X τ m−1 k+1 k+1 gik+1 − g̃ik+1 + gik − g̃ik εk+1 h + R , Y . i 2 i=1 (39) We have  Y k+1 ,Y k  η1k+1 η1j  Rk+1 , Y k+1  ≤ 1 ≤ 2  Y k+1 2 + Y 2 1 η1k+1 + ηij ≤ 2  (r1 ak + r2 bk )h2 Y k+1 L2 2 + 2 2 k 2   , , L2 Rk+1 4(r1 ak + r2 bk )h2 2 2 . Thus, Y k+1 1 ≤ 2 − × +  Y 2 2 k+1 2 2 + kY k22  − (r1 + r2 ) ×  X 1 k−1 [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] 2 j=0  η1k−j τL Y 2 2 h + (η1k+1 )2 h + ∆Y k−j k+1 2 2 + τL 4  Y k+1 2 2 + Y 18 2 η1k+1 2 2 k 2 2 h + ∆x Y + ∆Y k+1   2 2  k+1 2 2 /h2 + Rk+1 , Y k+1    2 /h  ≤ 1 2  Y k+1 2 2 + kY k22   2 2 (1 + ak )r1 + (1 + bk )r2 − × /h2 η1k+1 h + ∆x Y k+1 2 2     X 1 k−1 k−j 2 2 k−j 2 /h [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × − η1 h + ∆Y 2 2 j=0    2 2 τL 3τ L Y k+1 + Yk 2 2 4 4 2 2 (r1 ak + r2 bk )h L2 k+1 2 + Y Rk+1 + 2 2 2 2 L 4(r1 ak + r2 bk )h   2 1 Y k+1 + kY k22 ≤ 2 2     r1 + r2 2 k+1 2 k+1 2 − /h × h + ∆x Y η1 2 2     X 1 k−1 k−j 2 2 k−j 2 − /h [r1 (aj+1 − aj ) + r2 (bj+1 − bj )] × η1 h + ∆Y 2 2 j=0 + 2 2 τL r1 ak + r2 bk  k+1 2 3τ L + Y k+1 + Yk − h + ∆x Y k+1 η1 2 2 4 4 2 2 (r1 ak + r2 bk )h2 L2 k+1 k+1 2 + + . Y R 2 2 L2 4(r1 ak + r2 bk )h2  Lemma 3: Suppose Yk 2 v u m−1 u X = th |ηik |2 , Y k i=1 then Y k 2 2 ≤L Y k 2 ∞ ∞ = max 1≤i≤m−1 ηik , L2 ≤ 2 h|η1k |2 + ∆x Y k 2h  Proof : The first inequality is apparent. For the second inequality, let |ηik0 | = max |ηik k, 1≤i≤m−1 ηik0 = η1k + iX 0 −1 ∆x ηjk , ηik0 = − m−1 X j=i0 j=1 Thus 2|ηik0 | ≤ |η1k | + m−1 X j=1 19 |∆x ηjk |. ∆x ηjk . 2 2  . 2 2  /h2 (40) Using the Cauchy-Schwartz inequality, we have  4|ηik0 |2 ≤ 2m |η1k |2 + m−1 X j=1  |∆x ηjk |2  ≤ i 2L h 2 2 h|v | + k∆vk 1 2 . h2 Therefore, Yk 2 ∞ ≤ L2 h|η1k |2 + ∆x Y k 2h2  The lemma is proved. 2 2  . Applying Lemma 3, we obtain Y k+1 2 2 ≤ Yk + 2 2 + k X (r1 aj + r2 bj ) j=0 + k X (r1 aj + r2 bj ) j=0 3τ d1 Y k+1 2 2 2 +   τL Yk 2 2 2 2 η1k+1−j h + ∆x Y 2 η1k−j h + ∆x Y k−j + k+1−j 2 2 2 2  /h2 L2 Rk+1 2(r1 ak + r2 bk ) 2 2  2 /h   . (41) Let ρk = Y k 2 2 k X + (r1 aj + r2 bj ) j=0  2 η1k−j h + ∆x Y k−j 2 2  2 /h  (42) and using Lemma 1 Rk 2  ≤ C ak τ α + b k τ β   τ + h2 , r1 = Bτ β Aτ α , r2 = , Γ(α + 1) Γ(β + 1) then    2 1 3 1 − τ L ρk+1 ≤ 1 + τ L ρk + C ′ ak τ α + bk τ β τ + h2 . 2 2    Therefore, we obtain ρk+1 ρk+1 ≤    1 + 12 τ L ′ α β 2 2 ρ + C a τ + b τ . τ + h ≤ k k k 1 − 23 τ L 1 + 12 τ L 1 − 32 τ L  !k+1  ρ0 +  k X  C ′ aj τ α + b j τ β j=0 20  2 (43)  τ + h2  . (44) Note that ρ0 = 0, We can therefore conclude that there exists a positive constant Ĉ, such that k  X ρk+1 ≤ Ĉ ak τ α + b k τ β j=0 and k X  τ + h2 2 , ak τ α = (k + 1)α τ α ≤ T α . j=0 Thus  ρk+1 ≤ Ĉ T α + T β Because y K+1 2 2  τ + h2 2 . ≤ ρk+1 , we have Y k+1 2 2  ≤ Ĉ T α + T β  τ + h2 2 . Consequently, the following theorem of convergence is obtained. 2 Theorem 2: Let u(x, t) ∈ G(Ω) be the solution of (7-9) and assuming τ < 3L . Then the fractional implicit difference method defined by (25) is convergent, and there exists a positive constant C > 0 such that Y k+1 5 2 ≤ C(τ + h2 ). Numerical Results In this Section we illustrate some of the theory through numerical simulations. Example 1.We consider the following modified anomalous subdiffusion equation with a nonlinear source term: ∂ 1−α ∂ 1−β ∂u(x, t) = A 1−α + B 1−β ∂t ∂t ∂t !" ∂ 2 u(x, t) + f (x, t) + g(u(x, t), x, t), (45) ∂x2 # where A = B = 0.5, f (u, t) = 0, Γ(2 + α) 2α t ] Γ(1 + 2α) Γ(2 + β) 2β +ex [(1 + β)tβ − t ], Γ(1 + 2β) g(u(x, t), x, t) = ex [(1 + α)tα − with boundary condition and initial conditions 21 (46) 6 5 u(x,t) 4 3 2 1 0 1 0.8 1 0.6 0.8 0.6 0.4 0.4 0.2 0.2 0 Time−t (0−1) 0 Space−x (0−1) Fig. 1. Numerical solution of problem (45)-(47), when α = 0.5, β = 0.2 u(x, 0) = 0, 0 ≤ x ≤ 1, u(0, t) = t1+α + t1+β , u(1, t) = et1+α + et1+β . (47) The exact solution of Equations (45)-(47) is u(x, t) = ex (t1+α + t1+β ). We take α = 0.5, β = 0.2, 0 ≤ t ≤ 1, 0 ≤ x ≤ 1. The simulation results with α = 0.5, β = 0.2, 0 ≤ t ≤ 1, 0 ≤ x ≤ 1 are shown in Figure 1. The system exhibits behaviors of the solution and its derivatives of order α = 0.5, β = 0.2. We can also see that the u(x, t) increases with time. The comparisons of the numerical solution and exact solution are shown in Figure 2 when t = 0.2, 1, respectively. From Figure 2, it can be seen that the numerical solution is in good agreement with the exact solution. We take τ = 0.01, h = 0.1. The errors between the numerical solution and exact solution are shown in Table 1, when t = 1. From Table 1, we can see that the errors satisfy the relation Error ≤ C (τ + h2 ). Example 2. We consider the following modified anomalous subdiffusion equation with a nonlinear source term: ∂u(x, t) ∂ 1−α ∂ 1−β = A 1−α + B 1−β ∂t ∂t ∂t !" ∂ 2 u(x, t) + f (x, t) + g(u, x, t), ∂x2 22 # (48) Table 1 The error, numerical solution and exact solution, when t = 1, τ = 0.01, h = 0.1 Space Numerical solution Exact solution Error Error τ +h2 0.10 2.212387 2.210342 0.002045 0.10225 0.20 2.446536 2.442806 0.003730 0.18650 0.30 2.704757 2.699718 0.005039 0.25195 0.40 2.989597 2.983649 0.005947 0.29735 0.50 3.303859 3.297442 0.006417 0.32085 0.60 3.650636 3.644238 0.006398 0.31990 0.70 4.033338 4.027505 0.005833 0.29165 0.80 4.455724 4.451082 0.004642 0.23210 0.90 4.921941 4.919207 0.002734 0.13670 6 exact numerical 5 t=1.0 u(x,t)) 4 3 2 1 t=0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 Space−x (0−1) 0.7 0.8 0.9 1 Fig. 2. Numerical solution of problem (45)-(47), when α = 0.5, β = 0.2 with initial and boundary conditions:    8x, u(x, 0) =   −4x2 + 0 ≤ x ≤ 0.5, 22 x 3 + 43 , 0.5x ≤ x ≤ 2, u(0, t) = u(2, t) = 0, A = B = 1.0, f (x, t) = ex , g(u, x, t) = µ(u − u2 /K), 23 (49) 5 4 u(x,t) 3 2 1 0 0 0.5 2 1 1.5 0 1 0.5 0.5 1 Time−t (0−1) 0 Space−x (0−1) Fig. 3. Numerical solution of (48),(49), when α = 0.5, β = 0.9, t = 0 ∼ 1, 2 u(x,t=0.1) 1.5 1 0.5 0 1 0.8 2 0.6 1.5 0.4 1 0.2 Order−β (0−1) 0.5 0 0 Space−x (0−1) Fig. 4. Numerical solution of (48),(49), when α = 0.5, 0 ≤ β ≤ 1, t = 0.1, where g(u, x, t) is a Fisher nonlinear source term [1]. Here, we take µ = 0.5, K = 1. Figure 3 shows solution behaviors when α = 0.5, β = 0.9, t = 0 ∼ 1, while Figure 4 shows the response of the diffusion system for different real numbers 0 ≤ β ≤ 1, α = 0.5 at t = 0.4 and for different x. Figures 3-4 show that the system exhibits sub-diffusive behaviors and that the solution continuously depends on the time fractional derivative. 24 6 Conclusions In this paper, a new implicit numerical method for a modified anomalous subdiffusion equation with a nonlinear source term in a bounded domain has been described and demonstrated. We prove that the inplicit numerical method is stable and convergent using a new energy method. The implicit numerical method and analytical technique provide computationally effective tools for simulating the behavior of the solution of the modified anomalous subdiffusion equation with a nonlinear source term. This method and analytical technique can also be extended to any fractional integro-differential equations and higher-dimensional problems. Acknowledgements: Authors wish to thank the referees for their many constructive comments and suggestions to improve the paper. References [1] B. Baeumer, M.K. Mark and M. Meerschaert, Fractional reaction-diffusion equation for species growth and dispersal, Journal of Mathematical Biology,(2007), in press. [2] E. Brown, E. Wu, W. Zipfel, W. Webb, Measurement of molecular diffusion in solution by multiphoton fluorescence photobleaching recovery, Biophys. J., 77 (1999) 2837-2849. [3] A.V. Chechkin, R. Gorenflo, I.M. Sokolov, V.Yu. Gonchar, Distributed order time fractional diffusion equation, Frac. Calc. Appl. Anal., 6 (3) (2003) 259-279. [4] Chang-Ming Chen , F. Liu, I. Turner, and V. Anh , Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comp. Phys., 227 (2007) 886-897. [5] C. Chen, F. Liu and V. Anh, A Fourier method and an extrapolation technique for Stokes’ first problem for a heated generalized second grade fluid with fractional derivative, J. Comp. Appl. Math., (2008), in press (doi: 10.1016/j.cam.2008.03.01). 25 [6] D. Del-Castillo-Negrete, B.A. Carreras, V.E. Lynch, Front dynamics in reactiondiffusion systems with levy flights: a fractional diffusion approach, Phys. Rev. Lett. 91 (2003) 018302. [7] T. Feder, I. Brust-Mascher, J. Slattery, B. Baird, W. Webb, Constrained diffusion or immobile fraction on cell surfaces: a new interpretation, Biophys. J., 70 (1996) 2767-2773. [8] R. Ghosh, Mobility and clustering of individual low density lipoprotein receptor molecules on the surface of human skin fibroblasts, Ph.D. Thesis, Cornell University, Ithaca, NY, (1991). [9] R. Ghosh, W. Webb, Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules, Biophys. J., 66 (1994) 1301-1318. [10] B. Henry, S. Wearne, Fractional reaction-diffusion, Physica A 276 (2000) 448455. [11] T.A.M. Langlands and B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comp. Phys., 205 (2005) 719-736. [12] T.A.M. Langlands, Solution of a modified fractional diffusion equation, Physica A, 367 (2006) 136-144. [13] F. Liu, V. Anh and I. Turner, Numerical Solution of the Space Fractional Fokker-Planck Equation. J. Comp. Appl. Math., 166 (2004) 209-219. [14] J.P. Roop, Computational aspects of FEM approximation of fractional advection dispersion equation on bounded domains in R2 , J. Comp. Appl. Math., 193(1) (2006) 243-268. [15] F. Liu, P. Zhuang, V. Anh, I. Turner and K. Burrage, Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation, J. Appl. Math. Comp., 191 (2007) 12-21. 26 [16] F. Mainardi, M. Raberto, R. Gorenflo, E. Scalas, Fractional calculus and continuous-time finance II: the waiting-time distribution, Physica A, 287 (2000) 468-481. [17] J. Masoliver, M. Montero, Continuous-time random-walk model for financial distributions, Phys. Rev. E, 67 (2003) 021112. [18] J. Masoliver, M. Montero, J. Perello and G. Weiss, The CTRW in finance: direct and inverse problems (2003), hhttp://arxiv.org/abs/ cond-mat/0308017i. [19] J. Masoliver, M. Montero, J. Perello and G.H. Weiss, The continuous time random walk formalism in financial markets, J. Econ. Behav. Org., in press, hhttp://www.ffn.ub.es/papers/CTRW-aix.pdf and RePEc:sce:cplx03:24i. [20] M. Meerschaert and C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations. J. Comp. and Appl. Math., 172 (2004) 6577. [21] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep., 339 (2000) 1-77. [22] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A, 37 (2004) R161-R208. [23] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high frequency financial data: an empirical study, Physica A, 314 (2002) 749-755. [24] M. Saxton, Anomalous diffusion due to obstacles: a Monte Carlo Study, Biophys. J. 66 (1994) 394-401. [25] M. Saxton, Anomalous subdiffusion in fluorescence photobleaching recovery: a Monte Carlo study, Biophys. J., 81 (2001) 2226-2240. [26] M. Saxton, Anomalous diffusion due to binding: a Monte Carlo study, Biophys. J., 70 (1996) 1250-1262. [27] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance, Physica A, 284 (2000) 376-384. 27 [28] E. Scalas, R. Gorenflo, F. Mainardi, M. Mantelli, M. Raberto, Anomalous waiting times in high-frequency financial data 2003, hhttp://arxiv.org/abs/cond-mat/0310305i. [29] E. Scalas, Five years of continuous-time random walks in econophysics, in: A. Namatame (Ed.), Proceedings of WEHIA 2004, Kyoto, 2004. [30] E. Scalas, R. Gorenflo, F. Mainardi, M.M. Meerschaert, Speculative option valuation and the fractional diffusion equation, in: J. Sabatier, J. Tenreiro Machado (Eds.), Proceedings of the IFAC Workshop on Fractional Differentiation and its Applications, Bordeaux, (2004) 234-238. [31] E. Scalas, The application of continuous-time random walks in finance and economics, Physica A, 362(2) (2006) 225-239. [32] K. Seki, M. Wojcik, M. Tachiya, Fractional reaction-diffusion equation, J. Chem. Phys., 119 (4) (2003) 2165-2170. [33] K. Seki, M. Wojcik, M. Tachiya, Recombination kinetics in subdiffusive media, J. Chem. Phys., 119 (14) (2003) 7525-7533. [34] K. Seki, M. Wojcik and Tachiya, Fractional reaction-diffusion equation. J. Chem. Phys., 119 (2003) 2165-2174. [35] E. Sheets, G. Lee, R. Simson, K. Jacobson, Transient confinement of a glycosyl phosphatidylinositol-anchored protein in the plasma membrane, Biochemistry, 36 (1997) 12449-12458. [36] R. Simson, B. Yang, S. Moore, P. Doherty, F. Walsh, K. Jacobson, Structural mosaicism on the submicron scale in the plasma membrane, Biophys. J., 74 (1998) 297-308. [37] J. Slattery, Lateral mobility of FcRI on rat basophilic leukaemia cells as measured by single particle tracking using a novel bright fluorescent probe, Ph.D. Thesis, Cornell University, Ithaca, NY, (1991). [38] P. Smith, I. Morrison, K. Wilson, N. Fernandez, R. Cherry, Anomalous diffusion of major histocompatability complex class I molecules on HeLa cells determined by single particle tracking, Biophys. J., 76 (1999) 3331-3344. 28 [39] I.M. Sokolov, A.V. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta. Phys. Pol. B, 35 (2004) 1323. [40] I. Sokolov, J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s brownian motion, Chaos, 15 (2005) 026103. [41] W. C. Tan, T. Masuoka, Stokes’ first problem for a second grade fluid in a porous half-space with heated boundary, Int. J. Non-Linear Mech., 40 (2005) 515-522. [42] Q. Yu, F. Liu, V. Anh and I. Turner, Solving linear and nonlinear spacetime fractional reaction-diffusion equations by Adomian decomposition method, International J. for Numer. Meth. In Eng., 74 (2008) 138-158. [43] S.B. Yuste and K. Lindenberg, Subdiffusion limited reactions. Chem. Phys., 284 (2002) 169-180. [44] S.B. Yuste, L. Acedo and K. Lindenberg, Reaction front in an A + B → C reaction-subdiffusion process, Phys. Rev., E, 69 (2004) 036126. [45] S.B. Yuste and L. Acedo, An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion equations, SIAM J. Numer. Anal., 42(5) (2005) 1862-1874. [46] G. Zaslavsky, Fractional kinetic equation for Hamiltonian chaos. Chaotic advection, tracer dynamics and turbulent dispersion, Physica D, 76 (1994) 110122. [47] P. Zhuang, F. Liu, V. Anh and I. Turner, New solution and analytical techniques of the implicit numerical methods for the anomalous sub-diffusion equation, SIAM J. on Numerical Analysis, 46(2) (2008) 1079-1095. 29