International Journal of Differential Equations
Volume 2010 (2010), Article ID 464321, 22 pages
doi:10.1155/2010/464321
Research Article

Stability and Convergence of an Effective Numerical Method for the Time-Space Fractional Fokker-Planck Equation with a Nonlinear Source Term

Qianqian Yang, Fawang Liu, and Ian Turner

School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia

Received 25 May 2009; Revised 20 August 2009; Accepted 28 September 2009

Academic Editor: Om Agrawal

Copyright © 2010 Qianqian Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Fractional Fokker-Planck equations (FFPEs) have gained much interest recently for describing transport dynamics in complex systems that are governed by anomalous diffusion and nonexponential relaxation patterns. However, effective numerical methods and analytic techniques for the FFPE are still in their embryonic state. In this paper, we consider a class of time-space fractional Fokker-Planck equations with a nonlinear source term (TSFFPENST), which involve the Caputo time fractional derivative (CTFD) of order 𝛼 ∈ (0, 1) and the symmetric Riesz space fractional derivative (RSFD) of order 𝜇 ∈ (1, 2]. Approximating the CTFD and RSFD using the L1-algorithm and shifted Grünwald method, respectively, a computationally effective numerical method is presented to solve the TSFFPE-NST. The stability and convergence of the proposed numerical method are investigated. Finally, numerical experiments are carried out to support the theoretical claims.

1. Introduction

The Fokker-Planck equation (FPE) has commonly been used to describe the Brownian motion of particles. Normal diffusion in an external force field is often modeled in terms of the following Fokker-Planck equation (FPE) [1]: 𝜕 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 𝜕 𝑡 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 1 + 𝐾 1 𝜕 2 𝜕 𝑥 2 𝑢 ( 𝑥 , 𝑡 ) , ( 1 . 1 ) where 𝑚 is the mass of the diffusing test particle, 𝜂 1 denotes the fraction constant characterising the interaction between the test particle and its embedding, and the force is related to the external potential through 𝐹 ( 𝑥 ) = 𝑑 𝑉 ( 𝑥 ) / 𝑑 𝑥 . The FPE (1.1) is well studied for a variety of potential types, and the respective results have found wide application. In many studies of diffusion processes where the diffusion takes place in a highly nonhomogeneous medium, the traditional FPE may not be adequate [2, 3]. The nonhomogeneities of the medium may alter the laws of Markov diffusion in a fundamental way. In particular, the corresponding probability density of the concentration field may have a heavier tail than the Gaussian density, and its correlation function may decay to zero at a much slower rate than the usual exponential rate of Markov diffusion, resulting in long-range dependence. This phenomenon is known as anomalous diffusion [4]. Fractional derivatives play a key role in modeling particle transport in anomalous diffusion including the space fractional Fokker-Planck (advection-dispersion) equation describing Lévy flights, the time fractional Fokker-Planck equation depicting traps, and the time-space fractional equation characterizing the competition between Lévy flights and traps [5, 6]. Different assumptions on this probability density function lead to a variety of time-space fractional Fokker-Planck equations (TSFFPEs).

TSFFPE has been successfully used for modeling relevant physical processes. When the fractional differential equation is used to describe the asymptotic behavior of continuous time random walks, its solution corresponds to the Lévy walks, generalizing the Brownian motion to the Lévy motion. The following space fractional Fokker-Planck equation has been considered [2, 3, 7]: 𝜕 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 = − 𝑣 𝜕 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑥 + 𝐾 𝜇 𝑐 + 𝑎 𝐷 𝜇 𝑥 𝑢 ( 𝑥 , 𝑡 ) + 𝑐 − 𝑥 𝐷 𝜇 𝑏 , 𝑢 ( 𝑥 , 𝑡 ) ( 1 . 2 ) where 𝑣 is the drift of the process, that is, the mean advective velocity; 𝐾 𝜇 is the coefficient of dispersion; 𝑎 𝐷 𝜇 𝑥 and 𝑥 𝐷 𝜇 𝑏 are the left and right Riemann-Liouville space fractional derivatives of order 𝜇 given by 𝑎 𝐷 𝜇 𝑥 1 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 Γ ( 2 − 𝜇 ) 2 𝜕 𝑥 2 𝑥 𝑎 𝑢 ( 𝜉 , 𝑡 ) 𝑑 𝜉 ( 𝑥 − 𝜉 ) 𝜇 − 1 , 𝑥 𝐷 𝜇 𝑏 1 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 Γ ( 2 − 𝜇 ) 2 𝜕 𝑥 2 𝑏 𝑥 𝑢 ( 𝜉 , 𝑡 ) 𝑑 𝜉 ( 𝜉 − 𝑥 ) 𝜇 − 1 ; ( 1 . 3 ) 𝑐 + and 𝑐 − indicate the relative weight of transition probability; Benson et al. [2, 3] took 𝑐 + = 1 / 2 + 𝛽 / 2 and 𝑐 − = 1 / 2 − 𝛽 / 2 , ( − 1 ≤ 𝛽 ≤ 1 ) , which indicate the relative weight forward versus backward transition probability. If 𝑐 + = 𝑐 − = − 𝑐 𝜇 = − 1 / 2 c o s ( 𝜋 𝜇 / 2 ) , (1.2) can be rewritten in the following form: 𝜕 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 = − 𝑣 𝜕 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑥 + 𝐾 𝜇 𝜕 𝜇 𝑢 ( 𝑥 , 𝑡 ) 𝜕 | 𝑥 | 𝜇 , ( 1 . 4 ) where 𝜕 𝜇 / 𝜕 | 𝑥 | 𝜇 is the symmetric space fractional derivative of order 𝜇 ( 1 < 𝜇 ≤ 2 ). This is also referred to as the Riesz derivative [8], which contains a left Riemann-Liouville derivative ( 𝑎 𝐷 𝜇 𝑥 ) and a right Riemann-Liouville derivative ( 𝑥 𝐷 𝜇 𝑏 ), namely, 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) = − 𝑐 𝜇 𝑎 𝐷 𝜇 𝑥 + 𝑥 𝐷 𝜇 𝑏 𝑢 ( 𝑥 , 𝑡 ) . ( 1 . 5 )

As a model for subdiffusion in the presence of an external field, a time fractional extension of the FPE has been introduced as the time fractional Fokker-Planck equation (TFFPE) [5, 9]: 𝜕 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 𝑡 0 𝐷 𝑡 1 − 𝛼 𝜕 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝛼 𝜕 2 𝜕 𝑥 2 𝑢 ( 𝑥 , 𝑡 ) , ( 1 . 6 ) where the Riemann-Liouville operator 0 𝐷 𝑡 1 − 𝛼 , ( 0 < 𝛼 < 1 ) is defined through its operation: 0 𝐷 𝑡 1 − 𝛼 1 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 Γ ( 𝛼 ) 𝜕 𝑡 𝑡 0 𝑢 ( 𝑥 , 𝜂 ) ( 𝑡 − 𝜂 ) 1 − 𝛼 𝑑 𝜂 . ( 1 . 7 )

Yuste and Acedo [10] proposed an explicit finite difference method and a new von Neumann-type stability analysis for the anomalous subdiffusion equation (1.6) with 𝑉 ( 𝑥 ) = 0 . 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 L1-approximation scheme 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 et al. [12] presented a Fourier method for the anomalous subdiffusion equation, and they gave the stability analysis and the global accuracy analysis of the difference approximation scheme. Zhuang et al. [13] also proposed an implicit numerical method and an analytical technique for the anomalous subdiffusion equation. Chen et al. [14] proposed implicit and explicit numerical approximation schemes for the Stokes' first problem for a heated generalized second grade fluid with fractional derivatives. The stability and convergence of the numerical scheme 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 time-space fractional Fokker-Planck equation with a nonlinear source term are still in their infancy and are open problems.

Equation (1.6) can be written as the following equivalent form [15]: 0 𝐷 𝛼 𝑡 𝑢 ( 𝑥 , 𝑡 ) − 𝑢 ( 𝑥 , 0 ) 𝑡 − 𝛼 = 𝜕 Γ ( 1 − 𝛼 ) 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝛼 𝜕 2 𝜕 𝑥 2 𝑢 ( 𝑥 , 𝑡 ) . ( 1 . 8 ) By noting that [15] 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 0 𝐷 𝛼 𝑡 𝑢 ( 𝑥 , 𝑡 ) − 𝑢 ( 𝑥 , 0 ) 𝑡 − 𝛼 , Γ ( 1 − 𝛼 ) ( 1 . 9 ) we arrive at 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝛼 𝜕 2 𝜕 𝑥 2 𝑢 ( 𝑥 , 𝑡 ) , ( 1 . 1 0 ) where 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) / 𝜕 𝑡 𝛼 is the Caputo time fractional derivative (CTFD) of order 𝛼 ( 0 < 𝛼 < 1 ) with starting point at 𝑡 = 0 defined by [16] 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 1 Γ ( 1 − 𝛼 ) 𝑡 0 𝜕 𝑢 ( 𝑥 , 𝜂 ) 𝜕 𝜂 𝑑 𝜂 ( 𝑡 − 𝜂 ) 𝛼 . ( 1 . 1 1 ) The time-space fractional Fokker-Plank equation (TSFFPE), which describes the competition between subdiffusion and Lévy flights, is given by [5] 𝜕 𝑢 ( 𝑥 , 𝑡 ) = 𝜕 𝑡 0 𝐷 𝑡 1 − 𝛼 𝜕 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) , ( 1 . 1 2 ) or 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) , ( 1 . 1 3 ) where 𝐾 𝜇 𝛼 denotes the anomalous diffusion coefficient.

Schot et al. [17] investigated a fractional diffusion equation that employs time and space fractional derivatives by taking an absorbent (or source) term and an external force into account, which can be described by the following time-space fractional Fokker-Plank equation with an absorbent term and a linear external force: 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 𝜕 = − [ ] 𝜕 𝑥 𝐹 ( 𝑥 ) 𝑢 ( 𝑥 , 𝑡 ) + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) − 𝑡 0 𝑟 ( 𝑡 − 𝜂 ) 𝑢 ( 𝑥 , 𝜂 ) 𝑑 𝜂 , ( 1 . 1 4 ) where 𝐹 ( 𝑥 ) is the external force and 𝑟 ( 𝑡 ) is a time-dependent absorbent term, which may be related to a reaction diffusion process.

The fractional Fokker-Planck equations (FFPEs) have been recently treated by many authors and are presented as a useful approach for the description of transport dynamics in complex systems that are governed by anomalous diffusion and nonexponential relaxation patterns. The analytical solution of FFPE is only possible in simple and special cases [2, 3, 18] and the analytical solution provides a general representation in terms of Green's functions. We note that the representation of Green's functions is mostly expressed as convergent expansions in negative and positive power series. These special functions are not suitable for numerical evaluation when 𝑥 is sufficiently small or sufficiently large. Therefore, a new numerical strategy is important for solving these equations. Although numerical methods for the time fractional Fokker-Planck type equation, the space fractional Fokker-Plank type equation, and the time-space fractional Fokker-Planck type equation have been considered [7, 15, 19], numerical methods and stability and convergence analysis for the FFPE are quite limited and difficult. In fact, published papers on the numerical methods for the FFPE are sparse. We are unaware of any other published work on numerical methods for the time-space fractional Fokker-Planck type equation with a nonlinear source term. This motivates us to consider an effective numerical method for the time-space fractional Fokker-Planck equation with a nonlinear source term and to investigate its stability and convergence.

In this paper, we consider the following time-space fractional Fokker-Planck equation with a nonlinear source term (TSFFPE-NST): 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 𝑉 𝜕 𝑥 ( 𝑥 ) 𝑚 𝜂 𝛼 + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) + 𝑠 ( 𝑢 , 𝑥 , 𝑡 ) ( 1 . 1 5 ) subject to the boundary and initial conditions: 𝑢 𝑢 ( 𝑎 , 𝑡 ) = 𝑢 ( 𝑏 , 𝑡 ) = 0 , 0 ≤ 𝑡 ≤ 𝑇 , ( 𝑥 , 0 ) = 𝑢 0 ( 𝑥 ) , 𝑎 ≤ 𝑥 ≤ 𝑏 , ( 1 . 1 6 ) where 𝑝 ( 𝑥 ) = 𝑉 ( 𝑥 ) / 𝑚 𝜂 𝛼 is known as the drift coefficient. The nonlinear source (or absorbent) term 𝑠 ( 𝑢 , 𝑥 , 𝑡 ) is assumed to satisfy the Lipschitz condition: ‖ 𝑠 ( 𝑢 , 𝑥 , 𝑡 ) − 𝑠 ( 𝑣 , 𝑥 , 𝑡 ) ‖ ≤ 𝐿 ‖ 𝑢 − 𝑣 ‖ . ( 1 . 1 7 )

Let 𝑋 be a Banach space with associated norm ‖ 𝑢 ‖ . We say that 𝑠 ∶ 𝑋 → 𝑋 is globally Lipschitz continuous if for some 𝐿 > 0 , we have ‖ 𝑠 ( 𝑢 ) − 𝑠 ( 𝑣 ) ‖ ≤ 𝐿 ‖ 𝑢 − 𝑣 ‖ for all 𝑢 , 𝑣 ∈ 𝑋 , and is locally Lipschitz continuous, if the latter holds for ‖ 𝑢 ‖ , ‖ 𝑣 ‖ ≤ 𝑀 with 𝐿 = 𝐿 ( 𝑀 ) for any 𝑀 > 0 [20].

Let Ω = [ 𝑎 , 𝑏 ] × [ 0 , 𝑇 ] . In this paper, we suppose that the continuous problem (1.15)-(1.16) has a smooth solution 𝑢 ( 𝑥 , 𝑡 ) ∈ 𝐶 1 + 𝜇 , 2 𝑥 , 𝑡 ( Ω ) .

The rest of this paper is organized as follows. In Section 2, the Caputo time fractional derivative (CTFD) and the Riesz space fractional derivative (RSFD) are approximated by the L1-algorithm and the shifted Grünwald method, respectively. An effective numerical method (ENM) for solving the TSFFPE-NST (1.15)-(1.16) is proposed. The stability and convergence of the ENM are discussed in Sections 3 and 4, respectively. In Section 5, numerical experiments are carried out to support the theoretical analysis. Finally, some conclusions are drawn in Section 6.

2. An Effective Numerical Method for the TSFFPE-NST

In this section, we present an effective numerical method to simulate the solution behavior of the TSFFPE-NST (1.15)-(1.16). Let 𝑥 𝑙 = 𝑙 ℎ ( 𝑙 = 0 , 1 , … , 𝑀 ) and 𝑡 𝑛 = 𝑛 𝜏 ( 𝑛 = 0 , 1 , … , 𝑁 ), where ℎ = ( 𝑏 − 𝑎 ) / 𝑀 and 𝜏 = 𝑇 / 𝑁 are the spatial and temporal steps, respectively.

Firstly, adopting the L1-algorithm [21], we discretize the Caputo time fractional derivative as 𝜕 𝛼 𝑢 𝑥 , 𝑡 𝑛 + 1 𝜕 𝑡 𝛼 = 𝜏 − 𝛼 Γ ( 2 − 𝛼 ) 𝑛 𝑗 = 0 𝑏 𝑗 𝑢 𝑥 , 𝑡 𝑛 + 1 − 𝑗 − 𝑢 𝑥 , 𝑡 𝑛 − 𝑗 𝜏 + 𝑂 1 + 𝛼 , ( 2 . 1 ) where 𝑏 𝑗 = ( 𝑗 + 1 ) 1 − 𝛼 − 𝑗 1 − 𝛼 , 𝑗 = 0 , 1 , 2 , … , 𝑁 − 1 .

For the symmetric Riesz space fractional derivative, we use the following shifted Grünwald approximation [22]: 𝜕 𝜇 𝑢 𝑥 𝑙 , 𝑡 𝜕 | 𝑥 | 𝜇 ℎ = − − 𝜇 2 c o s ( 𝜋 𝜇 / 2 ) 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 − 𝑖 + 1 + , 𝑡 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 + 𝑖 − 1 ℎ , 𝑡 + 𝑂 𝑞 , ( 2 . 2 ) where the coefficients are defined by 𝑤 0 = 1 , 𝑤 𝑖 = ( − 1 ) 𝑖 𝜇 ( 𝜇 − 1 ) ⋯ ( 𝜇 − 𝑖 + 1 ) 𝑖 ! , 𝑖 = 1 , 2 , … , 𝑀 . ( 2 . 3 ) This formula is not unique because there are many different valid choices for 𝑤 𝑖 that lead to approximations of different orders 𝑞 [23]. The definition (2.2) provides order 𝑞 = 1 .

The first-order spatial derivative can be approximated by the backward difference scheme if 𝑝 ( 𝑥 ) < 0 , (otherwise, the forward difference scheme can be used if 𝑝 ( 𝑥 ) > 0 ): 𝜕 𝑝 𝑥 𝜕 𝑥 𝑙 𝑢 𝑥 𝑙 = 𝑝 𝑥 , 𝑡 𝑙 𝑢 𝑥 𝑙 𝑥 , 𝑡 − 𝑝 𝑙 − 1 𝑢 𝑥 𝑙 − 1 , 𝑡 ℎ + 𝑂 ( ℎ ) . ( 2 . 4 )

The nonlinear source term can be discretised either explicitly or implicitly. In this paper, we use an explicit method and evaluate the nonlinear source term at the previous time step: 𝑠 𝑢 𝑥 , 𝑡 𝑛 + 1 , 𝑥 , 𝑡 𝑛 + 1 𝑢 = 𝑠 𝑥 , 𝑡 𝑛 , 𝑥 , 𝑡 𝑛 + 𝑂 ( 𝜏 ) . ( 2 . 5 ) In this way, we avoid solving a nonlinear system at each time step and obtain an unconditionally stable and convergent numerical scheme, as shown in Section 3. However, the shortcoming of the explicit method is that it generates additional temporal error, as shown in (2.5).

Thus, using (2.1)–(2.5), we have 𝜏 − 𝛼 Γ ( 2 − 𝛼 ) 𝑛 𝑗 = 0 𝑏 𝑗 𝑢 𝑥 𝑙 , 𝑡 𝑛 + 1 − 𝑗 𝑥 − 𝑢 𝑙 , 𝑡 𝑛 − 𝑗 = 𝑝 𝑥 𝑙 𝑢 𝑥 𝑙 , 𝑡 𝑛 + 1 𝑥 − 𝑝 𝑙 − 1 𝑢 𝑥 𝑙 − 1 , 𝑡 𝑛 + 1 ℎ − 𝐾 𝜇 𝛼 ℎ − 𝜇 2 c o s ( 𝜋 𝜇 / 2 ) 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 − 𝑖 + 1 , 𝑡 𝑛 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 + 𝑖 − 1 , 𝑡 𝑛 + 1 𝑢 𝑥 + 𝑠 𝑙 , 𝑡 𝑛 , 𝑥 𝑙 , 𝑡 𝑛 𝜏 + 𝑂 1 + 𝛼 . + ℎ + 𝜏 ( 2 . 6 ) After some manipulation, (2.6) can be written in the following form: 𝑢 𝑥 𝑙 , 𝑡 𝑛 + 1 = 𝑏 𝑛 𝑢 𝑥 𝑙 , 𝑡 0 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑥 𝑙 , 𝑡 𝑛 − 𝑗 + 𝜇 0 ℎ 𝑝 𝑥 𝑙 𝑢 𝑥 𝑙 , 𝑡 𝑛 + 1 𝑥 − 𝑝 𝑙 − 1 𝑢 𝑥 𝑙 − 1 , 𝑡 𝑛 + 1 − 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 − 𝑖 + 1 , 𝑡 𝑛 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑥 𝑙 + 𝑖 − 1 , 𝑡 𝑛 + 1 + 𝜇 0 𝑠 𝑢 𝑥 𝑙 , 𝑡 𝑛 , 𝑥 𝑙 , 𝑡 𝑛 + 𝑅 𝑙 𝑛 + 1 , ( 2 . 7 ) where 𝜇 0 = 𝜏 𝛼 Γ ( 2 − 𝛼 ) > 0 , 𝑟 0 = 𝐾 𝜇 𝛼 ℎ − 𝜇 / 2 c o s ( 𝜋 𝜇 / 2 ) < 0 , and | | 𝑅 𝑙 𝑛 + 1 | | ≤ 𝐶 1 𝜏 𝛼 𝜏 1 + 𝛼 . + ℎ + 𝜏 ( 2 . 8 )

Let 𝑢 𝑛 𝑙 be the numerical approximation of 𝑢 ( 𝑥 𝑙 , 𝑡 𝑛 ) , and let 𝑠 𝑛 𝑙 be the numerical approximation of 𝑠 ( 𝑢 ( 𝑥 𝑙 , 𝑡 𝑛 ) , 𝑥 𝑙 , 𝑡 𝑛 ) . We obtain the following effective numerical method (ENM) of the TSFFPE-NST (1.15)-(1.16): 𝑢 𝑙 𝑛 + 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 + 𝜇 0 ℎ 𝑝 𝑙 𝑢 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝑢 𝑛 + 1 𝑙 − 1 − 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 + 𝑖 − 1 + 𝜇 0 𝑠 𝑛 𝑙 ( 2 . 9 ) for 𝑙 = 1 , 2 , … , 𝑀 − 1 , 𝑛 = 0 , 1 , 2 , … , 𝑁 − 1 . The boundary and initial conditions can be discretised using 𝑢 𝑛 0 = 𝑢 𝑛 𝑀 𝑢 = 0 , 𝑛 = 0 , 1 , 2 , … , 𝑁 , 0 𝑙 = 𝑢 0 ( 𝑙 ℎ ) , 𝑙 = 0 , 1 , 2 , … , 𝑀 . ( 2 . 1 0 )

Remark 2.1. If we use the implicit method to approximate the nonlinear source term, the numerical method of the TSFFPE-NST can be written as 𝑢 𝑙 𝑛 + 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 + 𝜇 0 ℎ 𝑝 𝑙 𝑢 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝑢 𝑛 + 1 𝑙 − 1 − 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 + 𝑖 − 1 + 𝜇 0 𝑠 𝑙 𝑛 + 1 , ( 2 . 1 1 ) that is, replace 𝑠 𝑛 𝑙 in (2.9) with 𝑠 𝑙 𝑛 + 1 . This numerical method is stable and convergent when the source term 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) satisfies the Lipschitz condition (1.17) (see, e.g., [24]).

Lemma 2.2 (see [19]). The coefficients 𝑏 𝑗 satisfy (1) 𝑏 𝑗 > 0 for 𝑗 = 0 , 1 , 2 , … , 𝑛 ; (2) 1 = 𝑏 0 > 𝑏 1 > ⋯ > 𝑏 𝑛 , 𝑏 𝑛 → 0 as 𝑛 → ∞ ; (3)when 0 < 𝛼 < 1 , l i m 𝑗 → ∞ 𝑏 𝑗 − 1 𝑗 𝛼 = l i m 𝑗 → ∞ 𝑗 − 1 1 + 𝑗 − 1 1 − 𝛼 = 1 − 1 . 1 − 𝛼 ( 2 . 1 2 ) Thus, there is a positive constant 𝐶 2 such that 𝑏 𝑗 − 1 ≤ 𝐶 2 𝑗 𝛼 , 𝑗 = 0 , 1 , 2 , … . ( 2 . 1 3 )

Lemma 2.3 (see [25]). The coefficients 𝑤 𝑖 satisfy (1) 𝑤 0 = 0 , 𝑤 1 = − 𝜇 < 0 , and 𝑤 𝑖 > 0 for 𝑖 = 2 , 3 , … , 𝑀 ; (2) ∑ ∞ 𝑖 = 0 𝑤 𝑖 = 0 , and ∑ 𝑛 𝑖 = 0 𝑤 𝑖 < 0 for ∀ 𝑛 ∈ ℕ .

3. Stability of the Effective Numerical Method

In this section, we analyze the stability of the ENM (2.9)-(2.10). Firstly, we rewrite (2.9) in the following form: 𝑢 𝑙 𝑛 + 1 − 𝜇 0 ℎ 𝑝 𝑙 𝑢 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝑢 𝑛 + 1 𝑙 − 1 + 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 + 𝑖 − 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 + 𝜇 0 𝑠 𝑛 𝑙 . ( 3 . 1 )

Let ̃ 𝑢 𝑛 𝑙 be the approximate solution of the ENM (3.1), and let ̃ 𝑠 𝑛 𝑙 be the approximation of 𝑠 𝑛 𝑙 . Setting 𝜌 𝑛 𝑙 = 𝑢 𝑛 𝑙 − ̃ 𝑢 𝑛 𝑙 , we obtain the following roundoff error equation: 𝜌 𝑙 𝑛 + 1 − 𝜇 0 ℎ 𝑝 𝑙 𝜌 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝜌 𝑛 + 1 𝑙 − 1 + 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝜌 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝜌 𝑛 + 1 𝑙 + 𝑖 − 1 = 𝑏 𝑛 𝜌 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝜌 𝑙 𝑛 − 𝑗 + 𝜇 0 𝑠 𝑛 𝑙 − ̃ 𝑠 𝑛 𝑙 ( 3 . 2 ) for 𝑙 = 1 , 2 , … , 𝑀 − 1 ; 𝑛 = 0 , 1 , … , 𝑁 − 1 .

We suppose that 𝑝 ( 𝑥 ) ≤ 0 and that 𝑝 ( 𝑥 ) decreases monotonically on [ 𝑎 , 𝑏 ] . This is based on the fact that physical considerations and stability dictate that 𝑝 ( 𝑥 ) < 0 [26, 27].

Assuming ‖ 𝜌 𝑛 ‖ ∞ = m a x 1 ≤ 𝑙 ≤ 𝑀 − 1 | 𝜌 𝑛 𝑙 | , and using mathematical induction, we obtain the following theorem.

Theorem 3.1. Suppose that 𝜌 𝑛 𝑙 ( 𝑙 = 1 , 2 , … , 𝑀 − 1 , 𝑛 = 1 , 2 , … , 𝑁 ) is the solution of the roundoff error equation (3.2), and the nonlinear source term 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) satisfies the Lipschitz condition (1.17), then there is a positive constant 𝐶 0 , such that ‖ 𝜌 𝑛 ‖ ∞ ≤ 𝐶 0 ‖ ‖ 𝜌 0 ‖ ‖ ∞ , 𝑛 = 1 , 2 , … , 𝑁 . ( 3 . 3 )

Proof. When 𝑛 = 1 , assume that | 𝜌 1 𝑙 0 | = m a x { | 𝜌 1 1 | , | 𝜌 1 2 | , … , | 𝜌 1 𝑀 − 1 | } . Because 𝑝 ( 𝑥 ) ≤ 0 and decreases monotonically on [ 𝑎 , 𝑏 ] , we have 𝜇 0 ≤ − 0 ℎ 𝑝 𝑙 0 − 𝑝 𝑙 0 − 1 | | | 𝜌 1 𝑙 0 | | | 𝜇 ≤ − 0 ℎ 𝑝 𝑙 0 | | | 𝜌 1 𝑙 0 | | | + 𝜇 0 ℎ 𝑝 𝑙 0 − 1 | | | 𝜌 1 𝑙 0 − 1 | | | . ( 3 . 4 ) Using the properties of 𝜔 𝑖 in Lemma 2.3, we have 0 ≤ 𝜇 0 𝑟 0 ⎡ ⎢ ⎢ ⎣ 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 | | | + 𝑀 − 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 | | | ⎤ ⎥ ⎥ ⎦ ≤ 𝜇 0 𝑟 0 ⎡ ⎢ ⎢ ⎣ 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 − 𝑖 + 1 | | | + 𝑀 − 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 + 𝑖 − 1 | | | ⎤ ⎥ ⎥ ⎦ . ( 3 . 5 ) Combining (3.4) with (3.5), using the Lipschitz condition (1.17) and smooth solution condition, we obtain | | | 𝜌 1 𝑙 0 | | | ≤ | | | 𝜌 1 𝑙 0 | | | − 𝜇 0 ℎ 𝑝 𝑙 0 | | | 𝜌 1 𝑙 0 | | | + 𝜇 0 ℎ 𝑝 𝑙 0 − 1 | | | 𝜌 1 𝑙 0 − 1 | | | + 𝜇 0 𝑟 0 ⎡ ⎢ ⎢ ⎣ 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 − 𝑖 + 1 | | | + 𝑀 − 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 | | | 𝜌 1 𝑙 0 + 𝑖 − 1 | | | ⎤ ⎥ ⎥ ⎦ ≤ | | | | | 𝜌 1 𝑙 0 − 𝜇 0 ℎ 𝑝 𝑙 0 𝜌 1 𝑙 0 + 𝜇 0 ℎ 𝑝 𝑙 0 − 1 𝜌 1 𝑙 0 − 1 + 𝜇 0 𝑟 0 ⎡ ⎢ ⎢ ⎣ 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 𝜌 1 𝑙 0 − 𝑖 + 1 + 𝑀 − 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 𝜌 1 𝑙 0 + 𝑖 − 1 ⎤ ⎥ ⎥ ⎦ | | | | | = | | | 𝑏 0 𝜌 0 𝑙 0 + 𝜇 0 𝑠 0 𝑙 0 − ̃ 𝑠 0 𝑙 0 | | | ≤ 𝑏 0 | | | 𝜌 0 𝑙 0 | | | + 𝜇 0 𝐿 | | | 𝜌 0 𝑙 0 | | | = 1 + 𝜇 0 𝐿 | | | 𝜌 0 𝑙 0 | | | . ( 3 . 6 ) Let 𝐶 = 1 + 𝜇 0 𝐿 . Thus, we obtain ‖ ‖ 𝜌 1 ‖ ‖ ∞ ‖ ‖ 𝜌 ≤ 𝐶 0 ‖ ‖ ∞ . ( 3 . 7 ) Now, suppose that ‖ ‖ 𝜌 𝑘 ‖ ‖ ∞ ‖ ‖ 𝜌 ≤ 𝐶 0 ‖ ‖ ∞ , 𝑘 = 2 , … , 𝑛 . ( 3 . 8 ) By assuming | 𝜌 𝑙 𝑛 + 1 0 | = m a x { | 𝜌 1 𝑛 + 1 | , | 𝜌 2 𝑛 + 1 | , … , | 𝜌 𝑛 + 1 𝑀 − 1 | } , we have that | | | 𝜌 𝑙 𝑛 + 1 0 | | | ≤ | | | | | 𝜌 𝑙 𝑛 + 1 0 − 𝜇 0 ℎ 𝑝 𝑙 0 𝜌 𝑙 𝑛 + 1 0 − 𝑝 𝑙 0 − 1 𝜌 𝑙 𝑛 + 1 0 − 1 + 𝜇 0 𝑟 0 ⎡ ⎢ ⎢ ⎣ 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 𝜌 𝑙 𝑛 + 1 0 − 𝑖 + 1 + 𝑀 − 𝑙 0 + 1 𝑖 = 0 𝑤 𝑖 𝜌 𝑙 𝑛 + 1 0 + 𝑖 − 1 ⎤ ⎥ ⎥ ⎦ | | | | | = | | | | | 𝑏 𝑛 𝜌 0 𝑙 0 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝜌 𝑙 𝑛 − 𝑗 0 + 𝜇 0 𝑠 𝑛 𝑙 0 − ̃ 𝑠 𝑛 𝑙 0 | | | | | ≤ 𝑏 𝑛 | | | 𝜌 0 𝑙 0 | | | + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 | | | 𝜌 𝑙 𝑛 − 𝑗 0 | | | + 𝜇 0 𝐿 | | | 𝜌 𝑛 𝑙 0 | | | . ( 3 . 9 ) Using (3.7) and (3.8), we have | | | 𝜌 𝑙 𝑛 + 1 0 | | | ≤ 𝑏 𝑛 | | | 𝜌 0 𝑙 0 | | | + 𝐶 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 | | | 𝜌 0 𝑙 0 | | | + 𝐶 𝜇 0 𝐿 | | | 𝜌 0 𝑙 0 | | | = 𝑏 𝑛 | | | 𝜌 0 𝑙 0 | | | 𝑏 + 𝐶 0 − 𝑏 𝑛 | | | 𝜌 0 𝑙 0 | | | + 𝐶 𝜇 0 𝐿 | | | 𝜌 0 𝑙 0 | | | = 𝑏 𝑛 𝜇 0 𝐿 + 𝐶 2 | | | 𝜌 0 𝑙 0 | | | . ( 3 . 1 0 ) Let 𝐶 0 = 𝑏 𝑛 𝜇 0 𝐿 + 𝐶 2 . Hence, we have ‖ ‖ 𝜌 𝑛 + 1 ‖ ‖ ∞ ≤ 𝐶 0 ‖ ‖ 𝜌 0 ‖ ‖ ∞ . ( 3 . 1 1 ) The proof of Theorem 3.1 is completed.

Applying Theorem 3.1, the following theorem of stability is obtained.

Theorem 3.2. Assuming that the nonlinear source term 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) satisfies the Lipschitz condition (1.17) and that the drift coefficient 𝑝 ( 𝑥 ) ≤ 0 decreases monotonically on [ 𝑎 , 𝑏 ] , the ENM defined by (2.9)-(2.10) is stable.

Remark 3.3. If 𝑝 ( 𝑥 ) > 0 and decreases monotonically on [ 𝑎 , 𝑏 ] , we can use the forward difference method to approximate the first-order spatial derivative and apply a similar analysis of stability.

Remark 3.4. In fact, for the case 𝑝 ( 𝑥 ) does not decrease monotonically, we can still obtain a stable numerical scheme by a minor change in our current ENM. We can expand the first term on the RHS of (1.15) as ( 𝜕 / 𝜕 𝑥 ) [ 𝑝 ( 𝑥 ) 𝑢 ( 𝑥 , 𝑡 ) ] = ( 𝑑 𝑝 / 𝑑 𝑥 ) 𝑢 ( 𝑥 , 𝑡 ) + 𝑝 ( 𝑥 ) ( 𝜕 𝑢 ( 𝑥 , 𝑡 ) / 𝜕 𝑥 ) , which enables us to group ( 𝑑 𝑝 / 𝑑 𝑥 ) 𝑢 ( 𝑥 , 𝑡 ) together with the nonlinear source term 𝑠 ( 𝑢 , 𝑥 , 𝑡 ) to obtain a new nonlinear source term 𝑠 ∗ ( 𝑢 , 𝑥 , 𝑡 ) = 𝑠 ( 𝑢 , 𝑥 , 𝑡 ) + ( 𝑑 𝑝 / 𝑑 𝑥 ) 𝑢 ( 𝑥 , 𝑡 ) . This way we can weaken the assumption on 𝑝 ( 𝑥 ) and the analysis given in this section still can be used.

Remark 3.5. If we use an implicit method to approximate the nonlinear source term, as shown in Remark 2.1, we can prove that the numerical method defined in (2.11) is stable when 1 − 𝜇 0 𝐿 > 0 , which is independent of the spatial step. In fact, when the time step is small, the condition 1 − 𝜇 0 𝐿 > 0 is generally satisfied.

4. Convergence of the Effective Numerical Method

In this section, we analyze the convergence of the ENM (2.9)-(2.10). Let 𝑢 ( 𝑥 𝑙 , 𝑡 𝑛 ) be the exact solution of the TSFFPE-NST (1.15)-(1.16) at mesh point ( 𝑥 𝑙 , 𝑡 𝑛 ) , and let 𝑢 𝑛 𝑙 be the numerical solution of the TSFFPE-NST (1.15)-(1.16) computed using the ENM (2.9)-(2.10). Define 𝜂 𝑛 𝑙 = 𝑢 ( 𝑥 𝑙 , 𝑡 𝑛 ) − 𝑢 𝑛 𝑙 and 𝐘 𝑛 = ( 𝜂 𝑛 1 , 𝜂 𝑛 1 , … , 𝜂 𝑛 𝑀 − 1 ) 𝑇 . Subtracting (2.9) from (2.7) leads to 𝜂 𝑙 𝑛 + 1 − 𝜇 0 ℎ 𝑝 𝑙 𝜂 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝜂 𝑛 + 1 𝑙 − 1 + 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝜂 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝜂 𝑛 + 1 𝑙 + 𝑖 − 1 = 𝑏 𝑛 𝜂 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝜂 𝑙 𝑛 − 𝑗 + 𝜇 0 𝑠 𝑢 𝑥 𝑙 , 𝑡 𝑛 , 𝑥 𝑙 , 𝑡 𝑛 − 𝑠 𝑛 𝑙 + 𝑅 𝑙 𝑛 + 1 , ( 4 . 1 ) where 𝑙 = 1 , 2 , … , 𝑀 − 1 ; 𝑛 = 0 , 1 , … , 𝑁 − 1 .

Assuming that ‖ 𝐘 𝑛 ‖ ∞ = m a x 1 ≤ 𝑙 ≤ 𝑀 − 1 | 𝜂 𝑛 𝑙 | and using mathematical induction, we obtain the following theorem.

Theorem 4.1. Assuming the nonlinear source term 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) satisfies the Lipschitz condition (1.17), and the drift coefficient 𝑝 ( 𝑥 ) ≤ 0 decreases monotonically on [ 𝑎 , 𝑏 ] , the ENM defined by (2.9)-(2.10) is convergent, and there exists a positive constant 𝐶 ∗ , such that ‖ 𝐘 𝑛 ‖ ∞ ≤ 𝐶 ∗ 𝜏 1 + 𝛼 + ℎ + 𝜏 , 𝑛 = 1 , 2 , … , 𝑁 . ( 4 . 2 )

Proof. Assume | 𝑅 𝑘 0 𝑙 0 | = m a x 1 ≤ 𝑙 ≤ 𝑀 − 1 , 1 ≤ 𝑛 ≤ 𝑁 | 𝑅 𝑛 𝑙 | . Following a similar argument to that presented above for the stability analysis of the ENM (2.9)-(2.10), when 𝑛 = 1 , assuming that | 𝜂 1 𝑙 0 | = m a x { | 𝜂 1 1 | , | 𝜂 1 2 | , … , | 𝜂 1 𝑀 − 1 | } , we have | | | 𝜂 1 𝑙 0 | | | ≤ | | | 𝑏 0 𝜂 0 𝑙 0 + 𝜇 0 𝑠 𝑢 𝑥 𝑙 0 , 𝑡 0 , 𝑥 𝑙 0 , 𝑡 0 − 𝑠 0 𝑙 0 + 𝑅 1 𝑙 0 | | | . ( 4 . 3 ) Utilising 𝐘 0 = 0 , the Lipschitz condition (1.17), and smooth solution condition, we obtain | | | 𝜂 1 𝑙 0 | | | ≤ 𝑏 0 | | | 𝜂 0 𝑙 0 | | | + 𝜇 0 𝐿 | | | 𝜂 0 𝑙 0 | | | + | | | 𝑅 𝑘 0 𝑙 0 | | | = | | | 𝑅 𝑘 0 𝑙 0 | | | . ( 4 . 4 ) Thus, ‖ ‖ 𝐘 1 ‖ ‖ ∞ ≤ | | | 𝑅 𝑘 0 𝑙 0 | | | . ( 4 . 5 ) Now, suppose that ‖ ‖ 𝐘 𝑘 ‖ ‖ ∞ ≤ 𝑏 − 1 𝑘 − 1 | | | 𝑅 𝑘 0 𝑙 0 | | | , 𝑘 = 1 , 2 , … , 𝑛 . ( 4 . 6 ) Using Lemma 2.2, 𝑏 𝑘 > 𝑏 𝑘 + 1 , we have ‖ ‖ 𝐘 𝑘 ‖ ‖ ∞ ≤ 𝑏 𝑛 − 1 | | | 𝑅 𝑘 0 𝑙 0 | | | . ( 4 . 7 ) Similarly, assuming | 𝜂 𝑙 𝑛 + 1 0 | = m a x { | 𝜂 1 𝑛 + 1 | , | 𝜂 2 𝑛 + 1 | , … , | 𝜂 𝑛 + 1 𝑀 − 1 | } , we have | | | 𝜂 𝑙 𝑛 + 1 0 | | | ≤ | | | | | 𝑏 𝑛 𝜂 0 𝑙 0 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝜂 𝑙 𝑛 − 𝑗 0 + 𝜇 0 𝑠 𝑢 𝑥 𝑙 0 , 𝑡 𝑛 , 𝑥 𝑙 0 , 𝑡 𝑛 − 𝑠 𝑛 𝑙 0 + 𝑅 𝑙 𝑛 + 1 0 | | | | | . ( 4 . 8 ) Utilising 𝐘 0 = 0 , the Lipschitz condition (1.17), and smooth solution condition, we obtain | | | 𝜂 𝑙 𝑛 + 1 0 | | | ≤ 𝑏 𝑛 − 1 𝑏 0 − 𝑏 𝑛 | | | 𝑅 𝑘 0 𝑙 0 | | | + 𝜇 0 𝐿 𝑏 𝑛 − 1 | | | 𝑅 𝑘 0 𝑙 0 | | | + | | | 𝑅 𝑘 0 𝑙 0 | | | = 𝑏 𝑛 − 1 𝑏 0 − 𝑏 𝑛 + 𝜇 0 𝐿 + 𝑏 𝑛 | | | 𝑅 𝑘 0 𝑙 0 | | | = 𝑏 𝑛 − 1 𝑏 0 + 𝜇 0 𝐿 | | | 𝑅 𝑘 0 𝑙 0 | | | = 𝐶 𝑏 𝑛 − 1 | | | 𝑅 𝑘 0 𝑙 0 | | | . ( 4 . 9 ) Hence, ‖ ‖ 𝐘 𝑛 + 1 ‖ ‖ ∞ ≤ 𝐶 𝑏 𝑛 − 1 | | | 𝑅 𝑘 0 𝑙 0 | | | . ( 4 . 1 0 ) Finally, utilising (2.8) and Lemma 2.2, 𝑏 𝑛 − 1 ≤ 𝐶 2 𝑛 𝛼 , we obtain the result on the convergence of the ENM (2.9)-(2.10), namely, ‖ 𝐘 𝑛 ‖ ∞ ≤ 𝐶 𝐶 1 𝐶 2 𝑛 𝛼 𝜏 𝛼 𝜏 1 + 𝛼 + ℎ + 𝜏 ≤ 𝐶 ∗ 𝜏 1 + 𝛼 + ℎ + 𝜏 ( 4 . 1 1 ) for 𝑛 = 1 , 2 , … , 𝑁 .

Remark 4.2. If we use an implicit method to approximate the nonlinear source term, as shown in Remark 2.1, we can prove that the numerical method defined in (2.11) is convergent when 1 − 𝜇 0 𝐿 > 0 , which is independent of the spatial step. In fact, when the time step is small, the condition 1 − 𝜇 0 𝐿 > 0 is generally satisfied.

5. Numerical Results

In this section, we present four numerical examples of the TSFFPE to demonstrate the accuracy of our theoretical analysis. We also use our solution method to illustrate the changes in solution behavior that arise when the exponent is varied from integer order to fractional order and to identify the differences between solutions with and without the external force term.

Example 5.1. Consider the following TSFFPE: 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 − 𝜐 𝜕 𝑥 + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 𝑢 𝑢 ( 𝑥 , 𝑡 ) + 𝑓 ( 𝑥 , 𝑡 ) , ( 𝑎 , 𝑡 ) = 𝑢 ( 𝑏 , 𝑡 ) = 0 , 0 ≤ 𝑡 ≤ 𝑇 , ( 𝑥 , 0 ) = 𝐾 𝜇 𝛼 ( 𝑥 − 𝑎 ) 2 ( 𝑏 − 𝑥 ) 2 , 𝑎 ≤ 𝑥 ≤ 𝑏 , ( 5 . 1 ) where 𝑓 ( 𝑥 , 𝑡 ) = ( 1 + 𝛼 ) 𝜐 Γ ( 1 + 𝛼 ) 𝑡 ( 𝑥 − 𝑎 ) 2 ( 𝑏 − 𝑥 ) 2 + 𝐾 𝜇 𝛼 𝐾 𝜇 𝛼 + 𝜐 𝑡 1 + 𝛼 [ 𝑔 ] 𝐾 2 c o s ( 𝜋 𝜇 / 2 ) ( 𝑥 − 𝑎 ) + 𝑔 ( 𝑏 − 𝑥 ) + 2 𝜐 𝜇 𝛼 + 𝜐 𝑡 1 + 𝛼 ( 𝑥 − 𝑎 ) ( 𝑏 − 𝑥 ) ( 𝑎 + 𝑏 − 2 𝑥 ) , 𝑔 ( 𝑥 ) = 4 ! Γ 𝑥 ( 5 − 𝜇 ) 4 − 𝜇 − 2 ( 𝑏 − 𝑎 ) 3 ! Γ 𝑥 ( 4 − 𝜇 ) 3 − 𝜇 + ( 𝑏 − 𝑎 ) 2 2 𝑥 Γ ( 3 − 𝜇 ) 2 − 𝜇 . ( 5 . 2 ) The exact solution of the TSFFPE (5.1) is found to be 𝐾 𝑢 ( 𝑥 , 𝑡 ) = 𝜇 𝛼 + 𝜐 𝑡 1 + 𝛼 ( 𝑥 − 𝑎 ) 2 ( 𝑏 − 𝑥 ) 2 , ( 5 . 3 ) which can be verified by direct fractional differentiation of the given solution, and substituting into the fractional differential equation.
In this example, we take 𝑎 = 0 , 𝑏 = 1 , 𝐾 𝜇 𝛼 = 2 5 , 𝜐 = 1 , 𝛼 = 0 . 8 , and 𝜇 = 1 . 9 . From Figure 1, it can be seen that the numerical solution using the ENM is in good agreement with the exact solution at different times 𝑇 , with ℎ = 1 / 4 0 and 𝜏 = 1 / 4 0 . The maximum errors of the ENM at time 𝑇 = 1 . 0 are presented in Table 1. It can be seen that the ENM is stable and convergent for solving the TSFFPE (5.1). The errors, as our theory indicated, satisfy the relationship e r r o r ≤ ( 𝜏 1 + 𝛼 + ℎ + 𝜏 ) .

tab1
Table 1: Maximum error behavior versus grid size reduction for Example 5.1 at 𝑇 = 1 . 0 .
464321.fig.001
Figure 1: Comparison of the numerical solution with the exact solution for Example 5.1 at 𝑇 = 1 . 0 , 3 . 0 , 5 . 0 . 𝑇 increases in the direction of the arrow.

Example 5.2. Consider the following TSFFPE-NST: 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝛾 𝑢 ( 𝑥 , 𝑡 ) − Γ ( 𝛽 ) 𝑡 0 ( 𝑡 − 𝜉 ) 𝛽 − 1 𝑢 ( 𝑥 , 𝜉 ) 𝑑 𝜉 , 𝑢 ( − 5 , 𝑡 ) = 𝑢 ( 5 , 𝑡 ) = 0 , 0 ≤ 𝑡 ≤ 𝑇 , 𝑢 ( 𝑥 , 0 ) = 𝛿 ( 𝑥 ) . ( 5 . 4 )
This example is a TSFFPE-NST without the external force term. In fact, it reduces to the fractional diffusion equation with an absorbent term. The formulae to approximate the absorbent term are presented in the appendix. Here, we take 𝛽 = 0 . 5 , 𝛾 = 1 , and 𝐾 𝜇 𝛼 = 1 . Figures 2–4 show the changes in the solution profiles of the TSFFPE-NST (5.4) when 𝛼 and 𝜇 are changed from integer to fraction at different times 𝑇 . We see that the solution profile of the fractional order model is characterized by a sharp peak and a heavy tail. The peak height in Figure 2 ( 𝛼 = 1 . 0 and 𝜇 = 2 . 0 ) decreases more rapidly than that in Figure 3 ( 𝛼 = 0 . 8 and 𝜇 = 1 . 8 ). Furthermore, when we choose 𝛼 = 0 . 5 and 𝜇 = 1 . 5 , a more interesting result can be observed; that is, the peak height in Figure 2 decreases more slowly than that shown in Figure 4 at the early time 𝑇 = 0 . 1 , but this trend reverses for the later times 𝑇 = 0 . 5 and 𝑇 = 1 . 0 . Hence, the TSFFPE-NST (5.4) may be useful to investigate several physical processes in the absence of an external force field by choosing appropriate 𝛼 and 𝜇 .

464321.fig.002
Figure 2: Numerical solutions for Example 5.2 with 𝛼 = 1 . 0 and 𝜇 = 2 . 0 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .
464321.fig.003
Figure 3: Numerical solutions for Example 5.2 with 𝛼 = 0 . 8 and 𝜇 = 1 . 8 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .
464321.fig.004
Figure 4: Numerical solutions for Example 5.2 with 𝛼 = 0 . 5 and 𝜇 = 1 . 5 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .

Example 5.3. Consider the following TSFFPE-NST: 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 𝜕 𝑥 𝑝 ( 𝑥 ) + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝛾 𝑢 ( 𝑥 , 𝑡 ) − Γ ( 𝛽 ) 𝑡 0 ( 𝑡 − 𝜉 ) 𝛽 − 1 𝑢 𝑢 ( 𝑥 , 𝜉 ) 𝑑 𝜉 , 𝑢 ( − 5 , 𝑡 ) = 𝑢 ( 5 , 𝑡 ) = 0 , 0 ≤ 𝑡 ≤ 𝑇 , ( 𝑥 , 0 ) = 𝛿 ( 𝑥 ) . ( 5 . 5 )
This example of the TSFFPE-NST incorporates the external force term with 𝑝 ( 𝑥 ) = − 1 and an absorbent term. The formula to approximate the absorbent term is presented in the appendix. Here, we take 𝛽 = 0 . 5 , 𝛾 = 1 , and 𝐾 𝜇 𝛼 = 1 . Figures 5–7 show the changes in the solution profiles of the TSFFPE-NST (5.5) when 𝛼 and 𝜇 are changed from integer order to fractional order at different times 𝑇 . Again, we see that the solution profile of the fractional order model is characterized by a sharp peak and a heavy tail. Furthermore, due to the presence of the external force term with 𝑝 ( 𝑥 ) = − 1 , the solution profiles are shifted to the right. It is worthwhile to note that the peak of the integer order model in Figure 5 ( 𝛼 = 1 . 0 and 𝜇 = 2 . 0 ) moves to the right as time increases, but the peak of the fractional order model in Figure 6 ( 𝛼 = 0 . 8 and 𝜇 = 1 . 8 ) and Figure 7 ( 𝛼 = 0 . 5 and 𝜇 = 1 . 5 ) does not move.
We also see that the peak heights in Figures 5 and 6 remain almost the same for increasing time. The peak height in Figure 5 decreases more slowly than that shown in Figure 7 at the early time 𝑇 = 0 . 1 , but this trend reverses for the later times 𝑇 = 0 . 5 and 𝑇 = 1 . 0 . Hence, the TSFFPE-NST (5.5) may be useful to investigate several physical processes within an external force field by choosing appropriate 𝛼 and 𝜇 .

464321.fig.005
Figure 5: Numerical solutions for Example 5.3 with 𝛼 = 1 . 0 and 𝜇 = 2 . 0 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .
464321.fig.006
Figure 6: Numerical solutions for Example 5.3 with 𝛼 = 0 . 8 and 𝜇 = 1 . 8 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .
464321.fig.007
Figure 7: Numerical solutions for Example 5.3 with 𝛼 = 0 . 5 and 𝜇 = 1 . 5 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .

Example 5.4. Consider the following TSFFPE-NST: 𝜕 𝛼 𝑢 ( 𝑥 , 𝑡 ) 𝜕 𝑡 𝛼 = 𝜕 𝜕 𝑥 𝑝 ( 𝑥 ) + 𝐾 𝜇 𝛼 𝜕 𝜇 𝜕 | 𝑥 | 𝜇 𝑢 ( 𝑥 , 𝑡 ) + 𝑟 𝑢 ( 𝑥 , 𝑡 ) 1 − 𝑢 ( 𝑥 , 𝑡 ) 𝐾 , 𝑢 ( 0 , 𝑡 ) = 𝑢 ( 5 , 𝑡 ) = 0 , 0 ≤ 𝑡 ≤ 𝑇 , 𝑢 ( 𝑥 , 0 ) = 𝑥 2 ( 5 − 𝑥 ) 2 , 0 ≤ 𝑥 ≤ 5 . ( 5 . 6 ) In applications to population biology, 𝑢 ( 𝑥 , 𝑡 ) is the population density at location 𝑥 ∈ ℛ and time 𝑡 > 0 . The nonlinear source term 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) = 𝑟 𝑢 ( 𝑥 , 𝑡 ) ( 1 − 𝑢 ( 𝑥 , 𝑡 ) / 𝐾 ) is Fisher's growth term that models population growth, where 𝑟 is the intrinsic growth rate of a species and 𝐾 is the environmental carrying capacity, representing the maximum sustainable population density [20, 28, 29].
In this example, we take 𝑟 = 0 . 2 , 𝐾 = 1 . Figure 8 shows the solution behavior when 𝛼 = 0 . 8 , 𝜇 = 1 . 6 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 , while Figure 9 shows the solution behavior with different values of 𝛼 between 0 and 1 and fixed value of 𝜇 = 1 . 8 at time 𝑇 = 1 . 0 . Figure 9 also shows that the system exhibits anomalous diffusion behavior and that the solution continuously depends on the time and space fractional derivatives. Although the source term for Fisher's equation 𝑠 ( 𝑢 ( 𝑥 , 𝑡 ) , 𝑥 , 𝑡 ) = 𝑟 𝑢 ( 𝑥 , 𝑡 ) ( 1 − 𝑢 ( 𝑥 , 𝑡 ) / 𝐾 ) is not globally Lipschitz continuous, the solution of the discrete numerical method still yields bounds on the solution of the continuous problem and the solution of the numerical method (ENM) converges to the unique solution of the continuous problem (5.6) as the time and space steps tend to zero [30].

464321.fig.008
Figure 8: Numerical solutions for Example 5.4 with 𝛼 = 0 . 8 and 𝜇 = 1 . 6 at different times 𝑇 = 0 . 1 , 0 . 5 , 1 . 0 .
464321.fig.009
Figure 9: Numerical solutions for Example 5.4 with fixed 𝜇 = 1 . 8 at time 𝑇 = 1 . 0 , and different values of 𝛼 = 0 . 3 , 0 . 4 , 0 . 5 , 0 . 6 , 0 . 7 , 0 . 8 , 0 . 9 , 1 . 0 .

6. Conclusions

In this paper, we have proposed an effective numerical method to solve the TSFFPE-NST and proved that the ENM is stable and convergent provided that the nonlinear source term satisfies the Lipschitz condition, the solution of the continuous problem satisfies the smooth solution condition, and 𝑝 ( 𝑥 ) can be either > 0 or < 0 . Numerical experiments have been carried out to support the theoretical claims. These numerical methods can also be used to investigate other types of fractional partial differential equations.

Appendix

Formulae for Examples 5.2 and 5.3

Let us start from (3.1), that is, 𝑢 𝑙 𝑛 + 1 − 𝜇 0 ℎ 𝑝 𝑙 𝑢 𝑙 𝑛 + 1 − 𝑝 𝑙 − 1 𝑢 𝑛 + 1 𝑙 − 1 + 𝜇 0 𝑟 0 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 + 𝑖 − 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 + 𝜇 0 𝑠 𝑛 𝑙 . ( A . 1 )

Now setting 𝑠 𝑛 𝑙 ∫ = − ( 𝛾 / Γ ( 𝛽 ) ) 𝑡 𝑛 0 ( 𝑡 𝑛 − 𝜉 ) 𝛽 − 1 𝑢 ( 𝑥 𝑙 , 𝜉 ) 𝑑 𝜉 , then we have 𝑠 𝑛 𝑙 𝛾 ≈ − Γ ( 𝛽 ) 𝑛 − 1 𝑗 = 0 𝑡 𝑗 + 1 𝑡 𝑗 𝑡 𝑛 − 𝜉 𝛽 − 1 𝑢 𝑥 𝑙 , 𝜉 𝑑 𝜉 . ( A . 2 )

Applying the Mean Value Theorem (M.V.T) for integration yields 𝑠 𝑛 𝑙 𝛾 ≈ − Γ ( 𝛽 ) 𝑛 − 1 𝑗 = 0 𝑢 𝑥 𝑙 , 𝜉 𝑗 𝑡 𝑗 + 1 𝑡 𝑗 𝑡 𝑛 − 𝜉 𝛽 − 1 𝑑 𝜉 , w h e r e 𝑡 𝑗 < 𝜉 𝑗 < 𝑡 𝑗 + 1 𝛾 ≈ − Γ ( 𝛽 ) 𝑛 − 1 𝑗 = 0 𝑢 𝑗 𝑙 + 𝑢 𝑙 𝑗 + 1 2 𝑡 𝑛 − 𝑡 𝑗 𝛽 𝛽 − 𝑡 𝑛 − 𝑡 𝑗 + 1 𝛽 𝛽 = − 𝛾 𝜏 𝛽 Γ ( 𝛽 ) ⋅ 2 ⋅ 𝛽 𝑛 − 1 𝑗 = 0 𝑢 𝑗 𝑙 + 𝑢 𝑙 𝑗 + 1 ( 𝑛 − 𝑗 ) 𝛽 − ( 𝑛 − 𝑗 − 1 ) 𝛽 = − 𝛾 𝜏 𝛽 2 Γ ( 𝛽 + 1 ) 𝑛 − 1 𝑗 = 0 𝑢 𝑙 𝑛 − 𝑗 − 1 + 𝑢 𝑙 𝑛 − 𝑗 ( 𝑗 + 1 ) 𝛽 − 𝑗 𝛽 = − 𝜇 1 𝑛 − 1 𝑗 = 0 𝑞 𝑗 𝑢 𝑙 𝑛 − 𝑗 − 1 + 𝑢 𝑙 𝑛 − 𝑗 , ( A . 3 ) where 𝜇 1 = 𝛾 𝜏 𝛽 / 2 Γ ( 𝛽 + 1 ) , 𝑞 𝑗 = ( 𝑗 + 1 ) 𝛽 − 𝑗 𝛽 , 𝑗 = 0 , 1 , … .

Also, we have 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 − 𝑖 + 1 + 𝑀 − 𝑙 + 1 𝑖 = 0 𝑤 𝑖 𝑢 𝑛 + 1 𝑙 + 𝑖 − 1 = 𝑀 − 1 𝑖 = 0 𝜂 𝑙 𝑖 𝑢 𝑖 𝑛 + 1 , ( A . 4 ) where 𝜂 𝑙 𝑖 = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑤 𝑙 − 𝑖 + 1 𝑤 , 1 ≤ 𝑖 ≤ 𝑙 − 2 , 0 + 𝑤 2 , 𝑖 = 𝑙 − 1 , 2 𝑤 1 𝑤 , 𝑖 = 𝑙 , 0 + 𝑤 2 𝑤 , 𝑖 = 𝑙 + 1 , 𝑖 − 𝑙 + 1 , 𝑙 + 2 ≤ 𝑖 ≤ 𝑀 − 1 . ( A . 5 )

Now, substituting (A.3) and (A.4) into (A.1), we obtain the numerical scheme for Example 5.2 as 𝑢 𝑙 𝑛 + 1 + 𝜇 0 𝑟 0 𝑀 − 1 𝑖 = 0 𝜂 𝑙 𝑖 𝑢 𝑖 𝑛 + 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 − 𝜇 1 𝑞 𝑗 𝑢 𝑙 𝑛 − 𝑗 − 1 + 𝑢 𝑙 𝑛 − 𝑗 , ( A . 6 ) and the numerical scheme for Example 5.3 as 𝜇 1 − 0 𝑝 𝑙 ℎ 𝑢 𝑙 𝑛 + 1 + 𝜇 0 𝑝 𝑙 − 1 ℎ 𝑢 𝑛 + 1 𝑙 − 1 + 𝜇 0 𝑟 0 𝑀 − 1 𝑖 = 0 𝜂 𝑙 𝑖 𝑢 𝑖 𝑛 + 1 = 𝑏 𝑛 𝑢 0 𝑙 + 𝑛 − 1 𝑗 = 0 𝑏 𝑗 − 𝑏 𝑗 + 1 𝑢 𝑙 𝑛 − 𝑗 − 𝜇 1 𝑞 𝑗 𝑢 𝑙 𝑛 − 𝑗 − 1 + 𝑢 𝑙 𝑛 − 𝑗 . ( A . 7 )

Acknowledgments

This research has been supported by a Ph.D. Fee Waiver Scholarship and a School of Mathematical Sciences Scholarship, QUT. The authors also wish to thank the referees for their constructive comments and suggestions.

References

  1. H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications, vol. 18 of Springer Series in Synergetics, Springer, Berlin, Germany, 2nd edition, 1989. View at MathSciNet
  2. D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert, “Application of a fractional advection-dispersion equation,” Water Resources Research, vol. 36, no. 6, pp. 1403–1412, 2000. View at Publisher · View at Google Scholar
  3. D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert, “The fractional-order governing equation of Lévy motion,” Water Resources Research, vol. 36, no. 6, pp. 1413–1423, 2000. View at Publisher · View at Google Scholar
  4. J.-P. Bouchaud and A. Georges, “Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications,” Physics Reports, vol. 195, no. 4-5, pp. 127–293, 1990. View at Publisher · View at Google Scholar · View at MathSciNet
  5. R. Metzler and J. Klafter, “Fractional Fokker-Planck equation: dispersive transport in an external force field,” Journal of Molecular Liquids, vol. 86, no. 1, pp. 219–228, 2000. View at Publisher · View at Google Scholar
  6. G. M. Zaslavsky, “Chaos, fractional kinetics, and anomalous transport,” Physics Reports, vol. 371, no. 6, pp. 461–580, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  7. F. Liu, V. Anh, and I. Turner, “Numerical solution of the space fractional Fokker-Planck equation,” Journal of Computational and Applied Mathematics, vol. 166, no. 1, pp. 209–219, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. R. Gorenflo and F. Mainardi, “Random walk models for space-fractional diffusion processes,” Fractional Calculus & Applied Analysis, vol. 1, no. 2, pp. 167–191, 1998. View at Zentralblatt MATH · View at MathSciNet
  9. C. W. Chow and K. L. Liu, “Fokker-Planck equation and subdiffusive fractional Fokker-Planck equation of bistable systems with sinks,” Physica A, vol. 341, no. 1–4, pp. 87–106, 2004. View at Publisher · View at Google Scholar · View at MathSciNet
  10. S. B. Yuste and L. Acedo, “An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations,” SIAM Journal on Numerical Analysis, vol. 42, no. 5, pp. 1862–1874, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  11. T. A. M. Langlands and B. I. Henry, “The accuracy and stability of an implicit solution method for the fractional diffusion equation,” Journal of Computational Physics, vol. 205, no. 2, pp. 719–736, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  12. C.-M. Chen, F. Liu, I. Turner, and V. Anh, “A Fourier method for the fractional diffusion equation describing sub-diffusion,” Journal of Computational Physics, vol. 227, no. 2, pp. 886–897, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  13. P. Zhuang, F. Liu, V. Anh, and I. Turner, “New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation,” SIAM Journal on Numerical Analysis, vol. 46, no. 2, pp. 1079–1095, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  14. C.-M. 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,” Journal of Computational and Applied Mathematics, vol. 223, no. 2, pp. 777–789, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  15. S. Chen, F. Liu, P. Zhuang, and V. Anh, “Finite difference approximations for the fractional Fokker-Planck equation,” Applied Mathematical Modelling, vol. 33, no. 1, pp. 256–273, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, vol. 198 of Mathematics in Science and Engineering, Academic Press, New York, NY, USA, 1999. View at MathSciNet
  17. A. Schot, M. K. Lenzi, L. R. Evangelista, L. C. Malacarne, R. S. Mendes, and E. K. Lenzi, “Fractional diffusion equation with an absorbent term and a linear external force: exact solution,” Physics Letters A, vol. 366, no. 4-5, pp. 346–350, 2007. View at Publisher · View at Google Scholar
  18. E. K. Lenzi, L. C. Malacarne, R. S. Mendes, and I. T. Pedron, “Anomalous diffusion, nonlinear fractional Fokker-Planck equation and solutions,” Physica A, vol. 319, pp. 245–252, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  19. P. Zhuang, F. Liu, V. Anh, and I. Turner, “Numerical treatment for the fractional Fokker-Planck equation,” The ANZIAM Journal, vol. 48, pp. C759–C774, 2006-2007. View at MathSciNet
  20. B. Baeumer, M. Kovács, and M. M. Meerschaert, “Numerical solutions for fractional reaction-diffusion equations,” Computers & Mathematics with Applications, vol. 55, no. 10, pp. 2212–2226, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  21. K. B. Oldham and J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, vol. 11 of Mathematics in Science and Engineering, Academic Press, New York, NY, USA, 1974. View at MathSciNet
  22. M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for two-sided space-fractional partial differential equations,” Applied Numerical Mathematics, vol. 56, no. 1, pp. 80–90, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  23. R. Lin and F. Liu, “Fractional high order methods for the nonlinear fractional ordinary differential equation,” Nonlinear Analysis: Theory, Methods & Applications, vol. 66, no. 4, pp. 856–869, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  24. F. Liu, C. Yang, and K. Burrage, “Numerical method and analytic technique of the modified anomalous subdiffusion equation with a nonlinear source term,” Journal of Computational and Applied Mathematics, vol. 231, no. 1, pp. 160–176, 2009. View at Publisher · View at Google Scholar · View at MathSciNet
  25. M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for fractional advection-dispersion flow equations,” Journal of Computational and Applied Mathematics, vol. 172, no. 1, pp. 65–77, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  26. W. Deng, “Finite element method for the space and time fractional Fokker-Planck equation,” SIAM Journal on Numerical Analysis, vol. 47, no. 1, pp. 204–226, 2008-2009. View at Publisher · View at Google Scholar · View at MathSciNet
  27. J. I. Ramos, “Damping characteristics of finite difference methods for one-dimensional reaction-diffusion equations,” Applied Mathematics and Computation, vol. 182, no. 1, pp. 607–609, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  28. D. del-Castillo-Negrete, B. A. Carreras, and V. E. Lynch, “Front dynamics in reaction-diffusion systems with Levy flights: a fractional diffusion approach,” Physical Review Letters, vol. 91, no. 1, Article ID 018302, 4 pages, 2003. View at Publisher · View at Google Scholar
  29. V. E. Lynch, B. A. Carreras, D. del-Castillo-Negrete, K. M. Ferreira-Mejias, and H. R. Hicks, “Numerical methods for the solution of partial differential equations of fractional order,” Journal of Computational Physics, vol. 192, no. 2, pp. 406–421, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  30. B. Baeumer, M. Kovács, and M. M. Meerschaert, “Fractional reproduction-dispersal equations and heavy tail dispersal kernels,” Bulletin of Mathematical Biology, vol. 69, no. 7, pp. 2281–2297, 2007. View at Publisher · View at Google Scholar · View at PubMed · View at MathSciNet