- Email: [email protected]
Method of approximate particular solutions for constant- and variable-order fractional diffusion models
Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Method of approximate particular solutions for constant- and variable-order fractional diffusion models Zhuo-Jia Fu a,b,n, Wen Chen a,nn, Leevan Ling b a b
College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China Department of Mathematics, Hong Kong Baptist University, Hong Kong, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 30 April 2014 Accepted 2 September 2014
The method of approximate particular solutions (MAPS) is an alternative radial basis function (RBF) meshless method, which is defined in terms of a linear combination of the particular solutions of the inhomogeneous governing equations with traditional RBFs as the source term. In this paper, we apply the MAPS to both constant- and variable-order time fractional diffusion models. In the discretization formulation, a finite difference scheme and the MAPS are used respectively to discretize time fractional derivative and spatial derivative terms. Numerical investigation examples show the present meshless scheme has highly accuracy and computationally efficiency for various fractional diffusion models. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Radial basis function Collocation method Meshless method Fractional diffusion
1. Introduction In recent decades, anomalous diffusion phenomena are extensively observed in a wide range of engineering and physics fields [1–4], such as contaminant transport, seepage, magnetic plasma, dissipation and turbulence. To describe anomalous diffusion phenomena, constant-order fractional diffusion equations are considered as recent alternative models and have received fantastic success [5–7]. However, various recent experimental results [8,9] show that constant-order fractional diffusion equations cannot fully capture some more complicated diffusion processes, whose diffusion behaviors depend on the time evolution, spatial variation or even concentration variation. To deal with these issues, variable-order fractional diffusion equations [10,11] have been introduced, in which the variable-order time fractional operator can be time-dependent, spatial-dependent, and/or concentration-dependent. Nowadays, finite difference methods (FDMs) are popular and dominant numerical techniques for temporal and spatial discretization of constant-order [12–16] and variable-order [17–20] fractional diffusion equations. Their convergence, accuracy, and stability have extensively been discussed in the literatures [21–24]. For the numerical simulations of constant-order fractional diffusion equations, with traditional FDMs for temporal discretization, several numerical methods have been introduced to spatial discretization of
n Corresponding author at: College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China. nn Corresponding author. E-mail addresses: [email protected] (Z.-J. Fu), [email protected] (W. Chen).
fractional derivative equations, such as the Fourier method [25], spectral method [26], finite element method [27–29], boundary element method [30], and radial basis function meshless collocation method [31–33]. In comparison with traditional FDMs for spatial discretization, these methods can reduce, to a certain extent, computing costs for large computational domain problems. In this work, we shall extend the idea to mitigate the computing costs in the numerical simulation of variable-order fractional diffusion equations. We will focus on constant- and variable-order time fractional diffusion equations, which only have fractional derivative in time and integer differential operator in space. We employ a finite difference method for temporal discretization and introduce an alternative radial basis function (RBF) meshless method, the method of approximate particular solutions (MAPS) [34–37], for spatial discretization. Chen et al. [38] first proposed the method of approximate particular solutions (MAPS) to solving partial differential equations. Then the MAPS has been successfully applied to various physical and engineering problems, such as anisotropic problems [39], nonlinear Poisson problems [40], wave problems [41], elasticity problems [42], Stokes flow problems [43], and convection-diffusion problems [44]. In comparison with the famous RBF method, also known as the Kansa method, the MAPS uses a newly derived RBF as interpolation basis function, which include some information from the considered governing equation operator. And some numerical experiments [45–47] demonstrate that the MAPS outperforms the Kansa method in terms of both the stability and accuracy, particularly in the evaluation of partial derivatives. This paper first applies the method of approximate particular solutions (MAPS), to 2D constant- and variable-order fractional diffusion problems. A brief outline of the paper is as follows. Section 2
http://dx.doi.org/10.1016/j.enganabound.2014.09.003 0955-7997/& 2014 Elsevier Ltd. All rights reserved.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
2
y
1.5
1
0.5
0
0
0.5
1 x
1.5
2
Fig. 1. Schematic configuration of uniform node distribution on a square domain (boundary nodes ‘o’ and inner nodes ‘n’).
Fig. 3. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(a)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
2. Methodology 2.1. Time fractional diffusion model Without loss of generality, we consider the following variableorder time fractional diffusion equations in a bounded domain Ω with piecewise smooth boundary ∂Ω ¼ Γ D þ Γ N (Γ D \ Γ N ¼ ∅) ∂αðtÞ uðx; t Þ ! ¼ DΔ þ v U ∇ λ uðx; t Þ þ Q ðx; t Þ; α ðt Þ ∂t 0 o αðt Þ o 1; x A Ω; t A ð0; T Þ;
ð1Þ
with boundary conditions Fig. 2. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(a)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
uðx; t Þ ¼ g 1 ðx; t Þ; ∂uðx; t Þ ¼ g 2 ðx; t Þ; ∂n
x A Γ D ; t A ð0; T Þ; x A Γ N ; t A ð0; T Þ;
ð2aÞ ð2bÞ
and initial condition describes the present computational formulations for fractional diffusion equations. In Section 3, the efficiency and accuracy of the present approach are examined with some benchmark examples. Finally, Section 4 concludes this paper with some remarks.
uðx; 0Þ ¼ u0 ðxÞ;
x A Ω;
ð3Þ
where Q ðx; t Þ, g 1 ðx; t Þ, g 2 ðx; t Þ and u0 ðxÞ are known functions; D the ! diffusion coefficient, λ the reaction coefficient, v the velocity vector, n the unit outward normal, T the total time to be considered,
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 4. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Fig. 5. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions. k
and ð∂αðtÞ =∂t αðtÞ Þ the variable-order time fractional derivative of order α(t) with respect to t defined by Z t ∂u x; η ∂αðt Þ uðx; t Þ 1 dη ¼ ð4Þ αðtÞ ; 0 o αðt Þ o 1: α ðt Þ ∂η ∂t 0 Γ 1α η t η For more details, readers are referred to Sun [11]. When α(t) is a constant and independent with t, the above-mentioned variableorder time fractional derivative goes back to the following famous Caputo constant-order time fractional derivative [48] Z t ∂u x; η ∂α uðx; t Þ 1 dη ¼ ð5Þ α ; 0 o α o 1 ∂t α ∂η Γ ð1 α Þ 0 t η Therefore, the constant-order time fractional derivative can be seen as a special case of the variable-order time fractional derivative. 2.2. Finite difference method for temporal discretization In the present method, we first introduce the finite difference method for temporal discretization of variable-order time fractional derivative (4). Then we have ∂αðt Þ uðx;t Þ ∂t αðt Þ
¼
Rt
∂uðx;ηÞ dη 1 0 Γ ð1 αðηÞÞ ∂η ðt ηÞαðt Þ
3
Z
j¼0
¼
ðj þ 1Þτ
∂u x; ξ 1 ∂ξ Γ 1 α t j þ 1 t
dξ ξ αðtj þ 1 Þ kþ1 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; k Z1
∑
jτ
> > : a0 u1 u0 ;
j¼1
ð6Þ
k¼0
where a0 ¼
τ αðt1 Þ τ αðtj þ 1 Þ ; ; aj ¼ Γ ð2 α ðt 1 ÞÞ Γ 2 α tj þ 1
h i bj ¼ ðj þ 1Þ1 αðtj þ 1 Þ j1 αðtj þ 1 Þ ;
j ¼ 1; 2; …; k:
More details can be found in Appendix A. For a constant-order time fractional derivative (5), we can still use formulation (6) and set αðt i Þ ¼ α; i ¼ 1; 2; …; k; k þ 1. Then the time fractional diffusion Eq. (1) can be rewritten as ! DΔuðx; t Þ þ v U ∇uðx; t Þ λuðx; t Þ þ Q ðx; t Þ 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; j¼1 ¼ > > : a0 u1 u0 ;
k Z1
ð7Þ
k¼0
Later, we use the θ-method to further temporal discretization of Eq. (7) and obtain
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
! DθΔuk þ 1 þ θ v U∇uk þ 1 λθ þ a0 uk þ 1 8 k > > > θ Q k þ 1 1 θ D Δ uk þ ! v U∇uk þ Q k λuk a0 uk þ ∑ aj bj uk j þ 1 uk j ; < j¼1 ¼ > ! > 1 0 0 0 0 0 > a D θ Q 1 θ Δ u þ v U∇u þQ λ u u ; : 0
where θ A ½0; 1, and method.
θ ¼0 is explicit method, θ ¼ 1 is implicit
2.3. Method of approximate particular solutions for spatial discretization The method of approximate particular solutions (MAPS) will be used for spatial discretization of Eq. (8), its computational formulation can be represented as uðxi ; t n þ 1 Þ ¼ ∑ αnj þ 1 Φ ‖xi xj ‖2 ; M
j¼1
ð9Þ
where ‖xi xj ‖2 the Euclidean distance between the collocation nodes xi and RBF centers xj , uðxi ; t n þ 1 Þ is the approximate solution
Fig. 6. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
kZ1
ð8Þ
k¼0
at the node xi at the instant t n þ 1 , and αnj þ 1 unknown coefficients at the instant t nþ 1 , M the total number of the collocation nodes, and Φ ‖xi xj ‖2 denotes the derived radial basis functions (RBFs) from Eq. (8). Since Eq. (8) can be considered as the standard inhomogeneous convection-diffusion equation with the governing ! operator Δ þ ð1=DÞ v U ∇ ððλθ þ a0 Þ=DθÞ, the present MAPS uses the following derived RBFs as interpolation basis function
Φðr Þ ¼
1
p=2
∑
μ2 i ¼ 0
Δ i p ðpÞ!!2 r lnðr Þ p þ 2 K 0 μr ; μ2 μ
ð10Þ
which are derived from the modified Helmholtz equation Δ μ2 ΦðrÞ ¼ ϕðr Þ with the thin plate spline RBF ðϕðr Þ ¼ r p lnðr Þ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi λθ þ a0 =ðDθÞ, K0 is the zero order p ¼ 2; 4; 6; …Þ in which μ ¼ modified Bessel function of the second kind. Here we list the first three derived RBF formulations from Eq. (10). When p ¼2, namely,
Fig. 7. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
ϕðr Þ ¼ r 2 lnðrÞ, we have
8 r2 lnðrÞ 4 lnðrÞ þ 4 < μ2 μ4 μ44 K 0 μr ; Φðr Þ ¼ : 44 þ 4γ4 þ 44 ln μ ; 2 μ μ μ
r a 0;
ð11aÞ
r ¼ 0;
when p ¼4, namely, ϕðr Þ ¼ r 4 lnðr Þ, we have 8 4 2 ðr Þ þ 1Þ 64 lnðr Þ þ 96 64K 0 ðμr Þ > < r μln2 ðrÞ 8r ð2 ln μ6 ; μ4 μ6 Φðr Þ ¼ μ 96 64γ 64 > : μ6 þ μ6 þ μ6 ln 2 ; when p ¼6, namely, ϕðr Þ ¼ r lnðr Þ, we have 8 n r6 lnðrÞ 12r4 ð3 lnðrÞ þ 1Þ 96r2 ð6 lnðrÞ þ 5Þ > μ2 ; > μ4 μ6 > > < o 2304 lnðrÞ þ 4224 2304K 0 ðμr Þ Φðr Þ ¼ μ8 μ8 > > > 4224 2304γ 2304 μ > : 8 þ 8 þ 8 ln ; 2 μ μ μ
r a 0; r ¼ 0;
5
calculated by the relative root mean square errors (RMSE) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi,sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 NT 1 NT 2 ~ i ; TÞ uðxi ; TÞÞ2 ∑ ðuðx ∑ u ðxi ; TÞ; RMSE ¼ NT i ¼ 1 NT i ¼ 1
ð12aÞ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ,sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ i ; TÞ ∂uðxi ; TÞ 2 1 NT ∂uðx 1 NT ∂uðxi ; TÞ 2 ; RMSEx ¼ ∑ ∑ NT i ¼ 1 ∂x1 ∂x1 NT i ¼ 1 ∂x1 ð12bÞ
ð11bÞ
~ i ; TÞ represent respectively the analytical and where uðxi ; TÞ and uðx numerical results evaluated at ðxi ; T Þ, and NT is the total number of uniform-distributed test nodes in the computational domain. Unless otherwise specified, NT¼ 441, T¼ 2 and θ ¼1 in all the following numerical cases.
6
r a 0; ð11cÞ r ¼ 0;
where Euler constant γ ¼ 0:57721566490153286….
3. Numerical results
Example 1. Consider the 2D constant-order fractional diffusion equation in Ω1 ∂α uðx; t Þ ¼ Δuðx; t Þ þ Q ðx; t Þ; ∂t α
0 o α o1; x A Ω; t A ð0; T Þ;
ð13Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions:
In this section, the efficiency, accuracy and convergence of the present method are tested to constant- and variable-order fractional diffusion equations under a square domain Ω1 with side length 2,
namely, Ω1 ¼ ðx1 ; x2 Þj 0 rx1 ; x2 r 2 . The numerical accuracy is
(a) Full Dirichlet boundary conditions
Fig. 8. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Fig. 9. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
uðx; t Þ ¼ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ;
x A Γ D ¼ ∂Ω; t A ð0; T Þ; ð14Þ
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
(b) Mixed boundary conditions 2
uðx; t Þ ¼ t ½x1 ð2 x1 Þ þx2 ð2 x2 Þ;
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ;
ð15aÞ
∂uðx; t Þ ∂ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ;
ð15bÞ
where Q ðx; t Þ ¼ ð2t 2 α =Γ ð3 αÞÞ þ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ þ 4t 2 . The exact solution is 2
uðx; t Þ ¼ t ½x1 ð2 x1 Þ þx2 ð2 x2 Þ;
x A Ω; t A ð0; T Þ:
ð16Þ
First, we set α ¼ 0.7 and place uniform nodes in the computational domain Ω1 as shown in Fig. 1. Figs. 2–7 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(a))–(11(c)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. Generally speaking, the numerical accuracy of the present method increases with the increasing interpolation node number M. However, with a large time step (dt ¼ 0.02, 0.01), the numerical accuracy first enhances with an increase of M, and then any further increase of M would not gain much improvement in terms of accuracy. This may be caused by the error generated from temporal discretization. From Figs. 2–7, it can be observed that the numerical accuracy of Example 1(a) is better than those of Example 1(b) by using the present method
Fig. 10. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
with different basis functions (11(a))–(11(c)). The present method with RBF formulations (11(b)) and (11(c)) can perform better than RBF formulation (11(a)). And it cannot provide the correct numerical results with the partial derivative term ðð∂uðx; t Þ=∂nÞ or ð∂uðx; t Þ=∂x1 ÞÞ by using RBF formulation (11(a)) for 2D fractional diffusion equations. Therefore, we only adopt RBF formulations (11 (b)) and (11(c)) as the basis functions in the MAPS to solve all the following numerical cases.
Example 2. Consider the 2D variable-order fractional diffusion equation with αðt Þ ¼ 0:8 þ 0:2t=T in Ω1 ∂αðtÞ uðx; t Þ ¼ Δuðx; t Þ þ Q ðx; t Þ; ∂t αðtÞ
0 o αðt Þ o1; x A Ω; t A ð0; T Þ;
ð17Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions: (a) Full Dirichlet boundary π x conditions 1 uðx; t Þ ¼ t 2 sin 2 π x sin 22 ; x A Γ D ¼ ∂Ω; t A ð0; TÞ; (b) Mixed boundary conditions π x π x 1 2 sin ; uðx; t Þ ¼ t 2 sin 2 2
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ;
ð19aÞ
Fig. 11. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
∂uðx; t Þ ∂ t 2 sin ðπ x1 =2Þ sin ðπ x2 =2Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ; ð19bÞ where Q ðx; t Þ ¼ ð2t 2 αðt Þ =Γ ð3 αðt ÞÞÞ þ ðπ 2 t 2 =2Þ sin ðπ x1 =2Þ sin ðπ x2 =2Þ . The exact solution is π x π x 1 2 sin ; x A Ω; t A ð0; T Þ: ð20Þ uðx; t Þ ¼ t 2 sin 2 2 Similar to the previous example, uniform nodes are placed in the computational domain Ω1 as shown in Fig. 1(a). Figs. 8–11 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(b))–(11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. Generally speaking, it can provide the acceptable numerical results (the error is less than 1%) by using the present method with node density (dh ¼1/8, namely, M ¼293) and dt r 0:02. With appropriate time step (dt ¼0.01), the numerical accuracy of the present method increase with the increasing interpolation node number M. However, with inappropriate time step (dt ¼0.02, 0.005, 0.002), the numerical accuracy first enhances with an increase of M, and then further increase of M may not gain much improvement in accuracy. It reveals that the present numerical solutions of variable-order time fractional diffusion models are more sensitive to the time discretization formulation and depend on the function type of the fractional derivative order. Therefore, coupling MAPS with other higher order time discretization formulations [49,50] will be a sensible way to improve the numerical accuracy in the
Fig. 12. Convergence rate (RMSE and RMSEx) of the present method with dt ¼ 0.02 by using different RBF formulations (11(b)) and (11(c)) in Example 3. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
7
solution of variable-order time fractional diffusion models. In this study, we will not further pursue this issue. Example 3. Consider the 2D variable-order fractional diffusion equation with αðt Þ ¼ 0:8 þ 0:2 sin 0:5π t=T in Ω1 ∂αðtÞ uðx; t Þ ¼ Δuðx; t Þ þQ ðx; t Þ; ∂t αðt Þ
0 o αðt Þ o 1; x A Ω; t A ð0; T Þ; ð21Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions: (a) Full Dirichlet boundary conditions π x π x 1 2 uðx; t Þ ¼ t 2 sin sin ; x A Γ D ¼ ∂Ω; t A ð0; T Þm 2 2 ð22Þ (b) Mixed boundary conditions π x π x 1 2 sin ; uðx; t Þ ¼ t 2 sin 2 2 t A ð0; T Þ;
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; ð23aÞ
∂uðx; t Þ ∂ t 2 sin ðπ x1 =2Þ sin ðπ x2 =2Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ; ð23bÞ 2 αðt Þ 2 2 =Γ ð3 αðt ÞÞÞ þ ðπ t =2Þ sin ðπ x1 =2Þ where Q ðx; tÞ ¼ ð2t sin ðπ x2 =2Þ . The exact solution is
Fig. 13. Convergence rate (RMSE and RMSEx) of the present method with dt ¼ 0.01 by using different RBF formulations (11(b)) and (11(c)) in Example 3. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
8
Table 1 Numerical errors (RMSE) by using the RBF formulation (11(b)) with both uniform and irregular nodes. Node distribution
RMSE
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
1.36E 04 4.50E 04
3.57E 03 6.38E 03
3.61E 04 1.90E 04
1.68E 03 2.53E 03
5.99E 03 6.17E 03
8.95E 03 1.07E 02
Table 2 Numerical errors (RMSEx) by using the RBF formulation (11(b)) with both uniform and irregular nodes. Node distribution
RMSEx
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
4.25E 03 4.72E 03
9.16E 03 3.15E 02
1.68E 02 1.38E 02
1.96E 02 6.23E 02
1.79E 02 1.61E 02
2.33E 02 6.18E 02
Table 3 Numerical errors (RMSE) by using the RBF formulation (11(c)) with both uniform and irregular nodes. Node distribution
RMSE
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
1.77E 04 1.89E 04
2.58E 03 3.69E 03
1.73E 04 1.62E 04
5.61E 04 7.94E 04
6.16E 03 6.17E 03
1.08E 02 1.10E 02
Table 4 Numerical errors (RMSEx) by using the RBF formulation (11(c)) with both uniform and irregular nodes. Node distribution
RMSEx
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
2.47E 03 1.65E 03
6.73E 03 1.72E 02
7.28E 03 7.64E 03
8.67E 03 2.72E 02
1.31E 02 1.40E 02
1.61E 02 3.16E 02
different time steps (dt ¼0.02, 0.01) in Example 3. From Figs. 12 and 13, we can see that the present method with node density (dh ¼1/8, namely, M¼293) provides the acceptable numerical results (the error is less than 1%). And it has the similar conclusion with Examples 1 and 2 that the numerical accuracy of Example 3 (a) is better than those of Example 3(b) by using the present method with different basis functions (11(b))–(11(c)).
1 α =0.3 α =0.5 α =0.7 α =0.9
u(-0.5,0,t)
0.8 0.6
Next we place irregular nodes in the computational domain Ω1 as shown in Fig. 13 to study the meshless feature of RBF-type methods. Tables 1–4 list the present numerical errors (RMSE and RMSEx) with the RBF formulations (11(b))–(11(c)) by using both uniform and irregular nodes with node number M¼484 and time step dt¼0.01. From Tables 1–4, it can be seen that the present method works equally well with irregular nodes. And the present method with RBF formulations (11(c)) can perform better than RBF formulation (11 (b)), especially for the partial derivative term ð∂uðxi ; TÞÞ=∂x1 .
0.4 0.2 0
0
0.5
1 t
1.5
2
Fig. 14. Temporal evolution of uð 0:5; 0; t Þ with different time fractional derivative order α by using the present method with RBF formulations (11(c)).
uðx; t Þ ¼ t 2 sin
π x 1
2
sin
π x 2
2
;
x A Ω; t A ð0; T Þ:
ð24Þ
Uniformly distributed nodes in Fig. 1(a) are used. Figs. 12 and 13 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(b))–(11(c)) by using
Example 4. The last case is the following fractional diffusion model in Ω1 ∂α uðx; t Þ ¼ Δuðx; t Þ; ∂t α
0 o α o 1; x A Ω; t A ð0; T Þ;
ð25Þ
with initial condition uðx; 0Þ ¼ f ðxÞ; x A Ω and under the mixed boundary conditions
uðx; t Þ ¼ 0; x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ; ð26aÞ
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
∂uðx; t Þ ¼ 0; ∂n where ( f ðxÞ ¼
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ;
2x1 þ2;
1 r x1 r1=2
ð2 2x1 Þ=3;
1=2 rx1 r 1
ð26bÞ
:
Finally, we set dt ¼0.01 and place uniform nodes (M¼488) in the computational domain Ω1 as shown in Fig. 1. Fig. 14 plots the temporal evolution of uð 0:5; 0; t Þ with different time fractional derivative order α by using the present method with RBF formulations (11(c)). Numerical investigation shows that the concentration u decays very fast at the beginning of diffusion process and it becomes more and more slowly at late times. This anomalous diffusion phenomenon becomes more and more obviously with the smaller time fractional derivative order α. 4. Conclusions This paper applies the method of approximate particular solutions (MAPS) to constant- and variable-order time fractional diffusion models. Discretization of the time-fractional and spatial derivatives is respectively done by a finite difference scheme and the MAPS. The present numerical experiments verify that the proposed method is a competitive meshless collocation method for constant- and variable-order time fractional diffusion models, and it provides satisfactory solutions and converges to the exact solutions with the increasing interpolation node number under
9
Moreover, it should be mentioned that the present numerical solutions of variable-order time fractional diffusion models are more sensitive to the time discretization formulation and depend on the function type of the fractional derivative order. Therefore, coupling MAPS with other higher order time discretization formulations [49,50] will be a sensible way to solve variable-order time fractional diffusion models. We leave the detailed numerical simulations to our future studies.
Acknowledgments The work described in this paper was supported by the National Science Funds of China (Grant nos. 11302069 and 11372097), the Fundamental Research Funds for the Central Universities, Hohai University (Grant no. 2013B32814), the National Basic Research Program of China (973 Project no. 2010CB832702), the National Science Funds for Distinguished Young Scholars of China (Grant no. 11125208), the 111 Project (Grant no. B12032), a CERG Grant of the Hong Kong Research Grant Council, and a FRG Grant of Hong Kong Baptist University.
Appendix A. Time discretization of variable-order time fractional derivative term The time discretization formulation (6) of variable-order time fractional derivative term (4) is derived as follows:
Z t ∂u x; η ∂αðt Þ uðx; t Þ 1 dη ¼ αðtÞ ∂η Γ ð1 α ðt ÞÞ 0 ∂t αðtÞ t η Z ðj þ 1Þτ k ∂u x; ξ 1 dξ ∑ α ðt j þ 1 Þ ∂ξ Γ 1 α t j þ 1 t j ¼ 0 jτ kþ1 ξ Z ðj þ 1Þτ k u x; t 1 dξ j þ 1 u x; t j ¼ ∑ τ Γ 1 α t j þ 1 t k þ 1 ξαðtj þ 1 Þ jτ j¼0 Z ðk j þ 1Þτ k u x; t 1 dη j þ 1 u x; t j ¼ ∑ τ Γ 1 α t k j þ 1 ηαðtk j þ 1 Þ ðk jÞτ j¼0 Z ðj þ 1Þτ k u x; t 1 dη k j þ 1 u x; t k j ¼ ∑ τ Γ 1 α t j þ 1 ηαðtj þ 1 Þ jτ j¼0 Z ðj þ 1Þτ k u x; t 1 dη k j þ 1 u x; t k j ¼ ∑ τ Γ 1 α t j þ 1 jτ ηαðtj þ 1 Þ j¼0
i h τ αðtj þ 1 Þ u x; t k j þ 1 u x; t k j ðj þ 1Þ1 αðt j þ 1 Þ j1 αðt j þ 1 Þ ; j ¼ 0 Γ 2 α tj þ 1 8 i α ðt j þ 1 Þ k > h > > τ αðt1 Þ uk þ 1 uk þ ∑ τ u x; t k j þ 1 u x; t k j ðj þ 1Þ1 αðtj þ 1 Þ j1 αðt j þ 1 Þ ; < Γ ð2 αðt 1 ÞÞ j ¼ 1Γ 2 α t j þ 1 ¼ > > > τ αðt1 Þ u1 u0 ; k ¼ 0 : Γ ð2 αðt 1 ÞÞ 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; k Z 1 j¼1 ¼ > > : a0 u1 u0 ; k ¼ 0 k
¼ ∑
both uniform and irregular distributions. Numerical investigations show that the present method with RBF formulation (11(c)) performs better than that with RBF formulations (11(a)) and (11 (b)), and it cannot provide the correct numerical results with the partial derivative term ðð∂uðx; t Þ=∂nÞ or ð∂uðx; t Þ=∂x1 ÞÞ by using RBF formulation (11(a)) for 2D fractional diffusion equations.
kZ1
where a0 ¼
τ αðt1 Þ τ α ðt j þ 1 Þ ; ; a ¼ Γ ð2 α ðt 1 ÞÞ j Γ 2 α t j þ 1
h i bj ¼ ðj þ 1Þ1 αðtj þ 1 Þ j1 αðtj þ 1 Þ ;
j ¼ 1; 2; …; k
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
10
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
References [1] Hilfer R. Applications of fractional calculus in physics. Singapore: World Scientific; 2000. [2] Metzler R, Klafter J. The random walk's guide to anomalous diffusion: a fractional dynamics approach. Phys Rep 2000;339(1):1–77. [3] Wyss W. The fractional diffusion equation. J Math Phys 1986;27:2782–5. [4] Gorenflo R, Luchko Y, Mainardi F. Wright functions as scale-invariant solutions of the diffusion-wave equation. J Comput Appl Math 2000;118(1–2):175–91. [5] Yang Q, Turner I, Liu F, Ilis M. Novel numerical methods for solving the timespace fractional diffusion equation in two dimensions. SIAM J Sci Comput 2011;33(3):1159–80. [6] Gorenflo R, Mainardi F, Moretti D, Pagnini G, Paradisi P. Discrete random walk models for space-time fractional diffusion. Chem Phys 2007;284(1–2):521–41. [7] Fu Z-J, Chen W, Yang H-T. Boundary particle method for Laplace transformed time fractional diffusion equations. J Comput Phys 2013;235:52–66. [8] Chechkin AV, Gorenflo R, Sokolov IM. Fractional diffusion in inhomogeneous media. J Phys A: Math Gen 2005;38(42):L679. [9] Sun H, Chen W, Chen Y. Variable-order fractional differential operators in anomalous diffusion modeling. Phys A: Stat Mech Appl 2009;388(21):4586–92. [10] Lorenzo CF, Hartley TT. Variable order and distributed order fractional operators. Nonlinear Dyn 2002;29(1–4):57–98. [11] Sun H, Chen W, Wei H. A comparative study of constant-order and variableorder fractional models in characterizing memory property of systems. Eur Phys J Spec Top 2011;193(1):185–92. [12] Tadjeran C, Meerschaert MM. A second-order accurate numerical method for the two-dimensional fractional diffusion equation. J Comput Phys 2007;220(2):813–23. [13] Murio DA. Implicit finite difference approximation for time fractional diffusion equations. Comput Math Appl 2008;56(4):1138–45. [14] Wang H, Wang K, Sircar T. A direct O(Nlog2N) finite difference method for fractional diffusion equations. J Comput Phys 2010;229(21):8095–104. [15] Ren J, Sun Z, Zhao X. Compact difference scheme for the fractional subdiffusion equation with Neumann boundary conditions. J Comput Phys 2013;232(1):456–67. [16] Zhao X, Sun Z. A box-type scheme for fractional sub-diffusion equation with Neumann boundary conditions. J Comput Phys 2011;230(15):6061–74. [17] Sun H, Chen W, Li C, Chen Y. Finite difference schemes for variable-order time fractional diffusion equation. Int J Bifurc Chaos 2012;22(04):1250085. [18] Shen S, Liu F, Chen J, Turner I, Anh V. Numerical techniques for the variable order time fractional diffusion equation. Appl Math Comput 2012;218(22):10861–70. [19] Chen CM, Liu F, Anh V, Turner I. Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation. SIAM J Sci Comput 2010;32(4):1740–60. [20] Lin R, Liu F, Anh V, Turner I. Stability and convergence of a new explicit finitedifference approximation for the variable-order nonlinear fractional diffusion equation. Appl Math Comput 2009;212(2):435–45. [21] Langlands TAM, Henry BI. The accuracy and stability of an implicit solution method for the fractional diffusion equation. J Comput Phys 2005;205(2): 719–36. [22] Cui M. Compact finite difference method for the fractional diffusion equation. J Comput Phys 2009;228(20):7792–804. [23] Zhuang P, Liu F. Implicit difference approximation for the time fractional diffusion equation. J Appl Math Comput 2006;22(3):87–99. [24] Zhang Y, Sun Z, Zhao X. Compact alternating direction implicit scheme for the two-dimensional fractional diffusion-wave equation. SIAM J Numer Anal 2012;50(3):1535–55. [25] Chen C-M, Liu F, Turner I, Anh V. A Fourier method for the fractional diffusion equation describing sub-diffusion. J Comput Phys 2007;227(2):886–97. [26] Lin Y, Xu C. Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys 2007;225(2):1533–52.
[27] Jiang Y, Ma J. High-order finite element methods for time-fractional partial differential equations. J Comput Appl Math 2011;235(11):3285–90. [28] Li C, Zhao Z, Chen Y. Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion. Comput Math Appl 2011;62(3):855–75. [29] Fix GJ, Roof JP. Least squares finite-element solution of a fractional order twopoint boundary value problem. Comput Math Appl 2004;48(7–8):1017–33. [30] Katsikadelis JT. The BEM for numerical solution of partial fractional differential equations. Comput Math Appl 2011;62(3):891–901. [31] Brunner H, Ling L, Yamamoto M. Numerical simulations of 2D fractional subdiffusion problems. J Comput Phys 2010;229(18):6613–22. [32] Chen W, Ye L, Sun H. Fractional diffusion equations by the Kansa method. Comput Math Appl 2010;59(5):1614–20. [33] Liu Q, Gu Y, Zhuang P, Liu F, Nie Y. An implicit RBF meshless approach for time fractional diffusion equations. Comput Mech 2011;48(1):1–12. [34] Chen CS, Fan CM, Wen PH. The method of approximate particular solutions for solving elliptic problems with variable coefficients. Int J Comput Methods 2011;8(03):545–59. [35] Ling L, Chen CS, Kwok TO. Adaptive method of particular solution For solving 3d inhomogeneous elliptic equations. Int J Comput Methods 2010;07(03):499–511. [36] Chen W, Fu ZJ, Chen CS. Recent advances on radial basis function collocation methods. Berlin: Springer; 2013. [37] Tsai C-C, Chen CS, Hsu T-W. The method of particular solutions for solving axisymmetric polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2009;33(12):1396–402. [38] Chen CS, Fan CM, Wen PH. The method of approximate particular solutions for solving certain partial differential equations. Numer Methods Partial Diff Equ 2012;28(2):506–22. [39] Zhu H. The method of approximate particular solutions for solving anisotropic elliptic problems. Eng Anal Bound Elem 2014;40:123–7. [40] Reutskiy SY. Method of particular solutions for nonlinear Poisson-type equations in irregular domains. Eng Anal Bound Elem 2013;37(2):401–8. [41] Kuo LH, Gu MH, Young DL, Lin CY. Domain type kernel-based meshless methods for solving wave equations. CMC-Comput Mater Contin 2013;33(3): 213–28. [42] Bustamante CA, Power H, Florez WF, Hang CY. The global approximate particular solution meshless method for two-dimensional linear elasticity problems. Int J Comput Math 2013;90(5):978–93. [43] Bustamante CA, Power H, Sua YH, Florez WF. A global meshless collocation particular solution method (integrated radial basis function) for twodimensional Stokes flow problems. Appl Math Modell 2013;37(6):4538–47. [44] Jiang T, Li M, Chen CS. The method of particular solutions for solving inverse problems of a nonhomogeneous convection-diffusion equation with variable coefficients. Numer Heat Transf, Part A: Appl 2012;61(5):338–52. [45] Yao G, Siraj ul I, Sarler B. Assessment of global and local meshless methods based on collocation with radial basis functions for parabolic partial differential equations in three dimensions. Eng Anal Bound Elem 2012;36(11): 1640–8. [46] Reutskiy SY. Method of particular solutions for solving PDEs of the second and fourth orders with variable coefficients. Eng Anal Bound Elem 2013;37(10): 1305–10. [47] Bustamante CA, Power H, Florez WF. A global meshless collocation particular solution method for solving the two-dimensional Navier–Stokes system of equations. Comput Math Appl 2013;65(12):1939–55. [48] Podlubny I. Fractional differential equations. San Diego: Academic Press; 1999. [49] Cuesta E, Lubich C, Palencia C. Convolution quadrature time discretization of fractional diffusion-wave equations. Math Comput 2006;75(254):673–96. [50] Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes. AIAA paper; 1981-1259.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Method of approximate particular solutions for constant- and variable-order fractional diffusion models Zhuo-Jia Fu a,b,n, Wen Chen a,nn, Leevan Ling b a b
College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China Department of Mathematics, Hong Kong Baptist University, Hong Kong, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 30 April 2014 Accepted 2 September 2014
The method of approximate particular solutions (MAPS) is an alternative radial basis function (RBF) meshless method, which is defined in terms of a linear combination of the particular solutions of the inhomogeneous governing equations with traditional RBFs as the source term. In this paper, we apply the MAPS to both constant- and variable-order time fractional diffusion models. In the discretization formulation, a finite difference scheme and the MAPS are used respectively to discretize time fractional derivative and spatial derivative terms. Numerical investigation examples show the present meshless scheme has highly accuracy and computationally efficiency for various fractional diffusion models. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Radial basis function Collocation method Meshless method Fractional diffusion
1. Introduction In recent decades, anomalous diffusion phenomena are extensively observed in a wide range of engineering and physics fields [1–4], such as contaminant transport, seepage, magnetic plasma, dissipation and turbulence. To describe anomalous diffusion phenomena, constant-order fractional diffusion equations are considered as recent alternative models and have received fantastic success [5–7]. However, various recent experimental results [8,9] show that constant-order fractional diffusion equations cannot fully capture some more complicated diffusion processes, whose diffusion behaviors depend on the time evolution, spatial variation or even concentration variation. To deal with these issues, variable-order fractional diffusion equations [10,11] have been introduced, in which the variable-order time fractional operator can be time-dependent, spatial-dependent, and/or concentration-dependent. Nowadays, finite difference methods (FDMs) are popular and dominant numerical techniques for temporal and spatial discretization of constant-order [12–16] and variable-order [17–20] fractional diffusion equations. Their convergence, accuracy, and stability have extensively been discussed in the literatures [21–24]. For the numerical simulations of constant-order fractional diffusion equations, with traditional FDMs for temporal discretization, several numerical methods have been introduced to spatial discretization of
n Corresponding author at: College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China. nn Corresponding author. E-mail addresses: [email protected] (Z.-J. Fu), [email protected] (W. Chen).
fractional derivative equations, such as the Fourier method [25], spectral method [26], finite element method [27–29], boundary element method [30], and radial basis function meshless collocation method [31–33]. In comparison with traditional FDMs for spatial discretization, these methods can reduce, to a certain extent, computing costs for large computational domain problems. In this work, we shall extend the idea to mitigate the computing costs in the numerical simulation of variable-order fractional diffusion equations. We will focus on constant- and variable-order time fractional diffusion equations, which only have fractional derivative in time and integer differential operator in space. We employ a finite difference method for temporal discretization and introduce an alternative radial basis function (RBF) meshless method, the method of approximate particular solutions (MAPS) [34–37], for spatial discretization. Chen et al. [38] first proposed the method of approximate particular solutions (MAPS) to solving partial differential equations. Then the MAPS has been successfully applied to various physical and engineering problems, such as anisotropic problems [39], nonlinear Poisson problems [40], wave problems [41], elasticity problems [42], Stokes flow problems [43], and convection-diffusion problems [44]. In comparison with the famous RBF method, also known as the Kansa method, the MAPS uses a newly derived RBF as interpolation basis function, which include some information from the considered governing equation operator. And some numerical experiments [45–47] demonstrate that the MAPS outperforms the Kansa method in terms of both the stability and accuracy, particularly in the evaluation of partial derivatives. This paper first applies the method of approximate particular solutions (MAPS), to 2D constant- and variable-order fractional diffusion problems. A brief outline of the paper is as follows. Section 2
http://dx.doi.org/10.1016/j.enganabound.2014.09.003 0955-7997/& 2014 Elsevier Ltd. All rights reserved.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
2
y
1.5
1
0.5
0
0
0.5
1 x
1.5
2
Fig. 1. Schematic configuration of uniform node distribution on a square domain (boundary nodes ‘o’ and inner nodes ‘n’).
Fig. 3. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(a)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
2. Methodology 2.1. Time fractional diffusion model Without loss of generality, we consider the following variableorder time fractional diffusion equations in a bounded domain Ω with piecewise smooth boundary ∂Ω ¼ Γ D þ Γ N (Γ D \ Γ N ¼ ∅) ∂αðtÞ uðx; t Þ ! ¼ DΔ þ v U ∇ λ uðx; t Þ þ Q ðx; t Þ; α ðt Þ ∂t 0 o αðt Þ o 1; x A Ω; t A ð0; T Þ;
ð1Þ
with boundary conditions Fig. 2. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(a)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
uðx; t Þ ¼ g 1 ðx; t Þ; ∂uðx; t Þ ¼ g 2 ðx; t Þ; ∂n
x A Γ D ; t A ð0; T Þ; x A Γ N ; t A ð0; T Þ;
ð2aÞ ð2bÞ
and initial condition describes the present computational formulations for fractional diffusion equations. In Section 3, the efficiency and accuracy of the present approach are examined with some benchmark examples. Finally, Section 4 concludes this paper with some remarks.
uðx; 0Þ ¼ u0 ðxÞ;
x A Ω;
ð3Þ
where Q ðx; t Þ, g 1 ðx; t Þ, g 2 ðx; t Þ and u0 ðxÞ are known functions; D the ! diffusion coefficient, λ the reaction coefficient, v the velocity vector, n the unit outward normal, T the total time to be considered,
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 4. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Fig. 5. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions. k
and ð∂αðtÞ =∂t αðtÞ Þ the variable-order time fractional derivative of order α(t) with respect to t defined by Z t ∂u x; η ∂αðt Þ uðx; t Þ 1 dη ¼ ð4Þ αðtÞ ; 0 o αðt Þ o 1: α ðt Þ ∂η ∂t 0 Γ 1α η t η For more details, readers are referred to Sun [11]. When α(t) is a constant and independent with t, the above-mentioned variableorder time fractional derivative goes back to the following famous Caputo constant-order time fractional derivative [48] Z t ∂u x; η ∂α uðx; t Þ 1 dη ¼ ð5Þ α ; 0 o α o 1 ∂t α ∂η Γ ð1 α Þ 0 t η Therefore, the constant-order time fractional derivative can be seen as a special case of the variable-order time fractional derivative. 2.2. Finite difference method for temporal discretization In the present method, we first introduce the finite difference method for temporal discretization of variable-order time fractional derivative (4). Then we have ∂αðt Þ uðx;t Þ ∂t αðt Þ
¼
Rt
∂uðx;ηÞ dη 1 0 Γ ð1 αðηÞÞ ∂η ðt ηÞαðt Þ
3
Z
j¼0
¼
ðj þ 1Þτ
∂u x; ξ 1 ∂ξ Γ 1 α t j þ 1 t
dξ ξ αðtj þ 1 Þ kþ1 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; k Z1
∑
jτ
> > : a0 u1 u0 ;
j¼1
ð6Þ
k¼0
where a0 ¼
τ αðt1 Þ τ αðtj þ 1 Þ ; ; aj ¼ Γ ð2 α ðt 1 ÞÞ Γ 2 α tj þ 1
h i bj ¼ ðj þ 1Þ1 αðtj þ 1 Þ j1 αðtj þ 1 Þ ;
j ¼ 1; 2; …; k:
More details can be found in Appendix A. For a constant-order time fractional derivative (5), we can still use formulation (6) and set αðt i Þ ¼ α; i ¼ 1; 2; …; k; k þ 1. Then the time fractional diffusion Eq. (1) can be rewritten as ! DΔuðx; t Þ þ v U ∇uðx; t Þ λuðx; t Þ þ Q ðx; t Þ 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; j¼1 ¼ > > : a0 u1 u0 ;
k Z1
ð7Þ
k¼0
Later, we use the θ-method to further temporal discretization of Eq. (7) and obtain
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
! DθΔuk þ 1 þ θ v U∇uk þ 1 λθ þ a0 uk þ 1 8 k > > > θ Q k þ 1 1 θ D Δ uk þ ! v U∇uk þ Q k λuk a0 uk þ ∑ aj bj uk j þ 1 uk j ; < j¼1 ¼ > ! > 1 0 0 0 0 0 > a D θ Q 1 θ Δ u þ v U∇u þQ λ u u ; : 0
where θ A ½0; 1, and method.
θ ¼0 is explicit method, θ ¼ 1 is implicit
2.3. Method of approximate particular solutions for spatial discretization The method of approximate particular solutions (MAPS) will be used for spatial discretization of Eq. (8), its computational formulation can be represented as uðxi ; t n þ 1 Þ ¼ ∑ αnj þ 1 Φ ‖xi xj ‖2 ; M
j¼1
ð9Þ
where ‖xi xj ‖2 the Euclidean distance between the collocation nodes xi and RBF centers xj , uðxi ; t n þ 1 Þ is the approximate solution
Fig. 6. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
kZ1
ð8Þ
k¼0
at the node xi at the instant t n þ 1 , and αnj þ 1 unknown coefficients at the instant t nþ 1 , M the total number of the collocation nodes, and Φ ‖xi xj ‖2 denotes the derived radial basis functions (RBFs) from Eq. (8). Since Eq. (8) can be considered as the standard inhomogeneous convection-diffusion equation with the governing ! operator Δ þ ð1=DÞ v U ∇ ððλθ þ a0 Þ=DθÞ, the present MAPS uses the following derived RBFs as interpolation basis function
Φðr Þ ¼
1
p=2
∑
μ2 i ¼ 0
Δ i p ðpÞ!!2 r lnðr Þ p þ 2 K 0 μr ; μ2 μ
ð10Þ
which are derived from the modified Helmholtz equation Δ μ2 ΦðrÞ ¼ ϕðr Þ with the thin plate spline RBF ðϕðr Þ ¼ r p lnðr Þ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi λθ þ a0 =ðDθÞ, K0 is the zero order p ¼ 2; 4; 6; …Þ in which μ ¼ modified Bessel function of the second kind. Here we list the first three derived RBF formulations from Eq. (10). When p ¼2, namely,
Fig. 7. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 1. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
ϕðr Þ ¼ r 2 lnðrÞ, we have
8 r2 lnðrÞ 4 lnðrÞ þ 4 < μ2 μ4 μ44 K 0 μr ; Φðr Þ ¼ : 44 þ 4γ4 þ 44 ln μ ; 2 μ μ μ
r a 0;
ð11aÞ
r ¼ 0;
when p ¼4, namely, ϕðr Þ ¼ r 4 lnðr Þ, we have 8 4 2 ðr Þ þ 1Þ 64 lnðr Þ þ 96 64K 0 ðμr Þ > < r μln2 ðrÞ 8r ð2 ln μ6 ; μ4 μ6 Φðr Þ ¼ μ 96 64γ 64 > : μ6 þ μ6 þ μ6 ln 2 ; when p ¼6, namely, ϕðr Þ ¼ r lnðr Þ, we have 8 n r6 lnðrÞ 12r4 ð3 lnðrÞ þ 1Þ 96r2 ð6 lnðrÞ þ 5Þ > μ2 ; > μ4 μ6 > > < o 2304 lnðrÞ þ 4224 2304K 0 ðμr Þ Φðr Þ ¼ μ8 μ8 > > > 4224 2304γ 2304 μ > : 8 þ 8 þ 8 ln ; 2 μ μ μ
r a 0; r ¼ 0;
5
calculated by the relative root mean square errors (RMSE) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi,sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 NT 1 NT 2 ~ i ; TÞ uðxi ; TÞÞ2 ∑ ðuðx ∑ u ðxi ; TÞ; RMSE ¼ NT i ¼ 1 NT i ¼ 1
ð12aÞ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ,sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ i ; TÞ ∂uðxi ; TÞ 2 1 NT ∂uðx 1 NT ∂uðxi ; TÞ 2 ; RMSEx ¼ ∑ ∑ NT i ¼ 1 ∂x1 ∂x1 NT i ¼ 1 ∂x1 ð12bÞ
ð11bÞ
~ i ; TÞ represent respectively the analytical and where uðxi ; TÞ and uðx numerical results evaluated at ðxi ; T Þ, and NT is the total number of uniform-distributed test nodes in the computational domain. Unless otherwise specified, NT¼ 441, T¼ 2 and θ ¼1 in all the following numerical cases.
6
r a 0; ð11cÞ r ¼ 0;
where Euler constant γ ¼ 0:57721566490153286….
3. Numerical results
Example 1. Consider the 2D constant-order fractional diffusion equation in Ω1 ∂α uðx; t Þ ¼ Δuðx; t Þ þ Q ðx; t Þ; ∂t α
0 o α o1; x A Ω; t A ð0; T Þ;
ð13Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions:
In this section, the efficiency, accuracy and convergence of the present method are tested to constant- and variable-order fractional diffusion equations under a square domain Ω1 with side length 2,
namely, Ω1 ¼ ðx1 ; x2 Þj 0 rx1 ; x2 r 2 . The numerical accuracy is
(a) Full Dirichlet boundary conditions
Fig. 8. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Fig. 9. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(b)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
uðx; t Þ ¼ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ;
x A Γ D ¼ ∂Ω; t A ð0; T Þ; ð14Þ
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
(b) Mixed boundary conditions 2
uðx; t Þ ¼ t ½x1 ð2 x1 Þ þx2 ð2 x2 Þ;
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ;
ð15aÞ
∂uðx; t Þ ∂ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ;
ð15bÞ
where Q ðx; t Þ ¼ ð2t 2 α =Γ ð3 αÞÞ þ t 2 ½x1 ð2 x1 Þ þ x2 ð2 x2 Þ þ 4t 2 . The exact solution is 2
uðx; t Þ ¼ t ½x1 ð2 x1 Þ þx2 ð2 x2 Þ;
x A Ω; t A ð0; T Þ:
ð16Þ
First, we set α ¼ 0.7 and place uniform nodes in the computational domain Ω1 as shown in Fig. 1. Figs. 2–7 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(a))–(11(c)) by using different time steps (dt ¼0.02, 0.01, 0.005, 0.002) in Example 1. Generally speaking, the numerical accuracy of the present method increases with the increasing interpolation node number M. However, with a large time step (dt ¼ 0.02, 0.01), the numerical accuracy first enhances with an increase of M, and then any further increase of M would not gain much improvement in terms of accuracy. This may be caused by the error generated from temporal discretization. From Figs. 2–7, it can be observed that the numerical accuracy of Example 1(a) is better than those of Example 1(b) by using the present method
Fig. 10. Convergence rate (RMSE) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
with different basis functions (11(a))–(11(c)). The present method with RBF formulations (11(b)) and (11(c)) can perform better than RBF formulation (11(a)). And it cannot provide the correct numerical results with the partial derivative term ðð∂uðx; t Þ=∂nÞ or ð∂uðx; t Þ=∂x1 ÞÞ by using RBF formulation (11(a)) for 2D fractional diffusion equations. Therefore, we only adopt RBF formulations (11 (b)) and (11(c)) as the basis functions in the MAPS to solve all the following numerical cases.
Example 2. Consider the 2D variable-order fractional diffusion equation with αðt Þ ¼ 0:8 þ 0:2t=T in Ω1 ∂αðtÞ uðx; t Þ ¼ Δuðx; t Þ þ Q ðx; t Þ; ∂t αðtÞ
0 o αðt Þ o1; x A Ω; t A ð0; T Þ;
ð17Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions: (a) Full Dirichlet boundary π x conditions 1 uðx; t Þ ¼ t 2 sin 2 π x sin 22 ; x A Γ D ¼ ∂Ω; t A ð0; TÞ; (b) Mixed boundary conditions π x π x 1 2 sin ; uðx; t Þ ¼ t 2 sin 2 2
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ;
ð19aÞ
Fig. 11. Convergence rate (RMSEx) of the present method with the derived RBF formulation (11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
∂uðx; t Þ ∂ t 2 sin ðπ x1 =2Þ sin ðπ x2 =2Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ; ð19bÞ where Q ðx; t Þ ¼ ð2t 2 αðt Þ =Γ ð3 αðt ÞÞÞ þ ðπ 2 t 2 =2Þ sin ðπ x1 =2Þ sin ðπ x2 =2Þ . The exact solution is π x π x 1 2 sin ; x A Ω; t A ð0; T Þ: ð20Þ uðx; t Þ ¼ t 2 sin 2 2 Similar to the previous example, uniform nodes are placed in the computational domain Ω1 as shown in Fig. 1(a). Figs. 8–11 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(b))–(11(c)) by using different time steps (dt ¼ 0.02, 0.01, 0.005, 0.002) in Example 2. Generally speaking, it can provide the acceptable numerical results (the error is less than 1%) by using the present method with node density (dh ¼1/8, namely, M ¼293) and dt r 0:02. With appropriate time step (dt ¼0.01), the numerical accuracy of the present method increase with the increasing interpolation node number M. However, with inappropriate time step (dt ¼0.02, 0.005, 0.002), the numerical accuracy first enhances with an increase of M, and then further increase of M may not gain much improvement in accuracy. It reveals that the present numerical solutions of variable-order time fractional diffusion models are more sensitive to the time discretization formulation and depend on the function type of the fractional derivative order. Therefore, coupling MAPS with other higher order time discretization formulations [49,50] will be a sensible way to improve the numerical accuracy in the
Fig. 12. Convergence rate (RMSE and RMSEx) of the present method with dt ¼ 0.02 by using different RBF formulations (11(b)) and (11(c)) in Example 3. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
7
solution of variable-order time fractional diffusion models. In this study, we will not further pursue this issue. Example 3. Consider the 2D variable-order fractional diffusion equation with αðt Þ ¼ 0:8 þ 0:2 sin 0:5π t=T in Ω1 ∂αðtÞ uðx; t Þ ¼ Δuðx; t Þ þQ ðx; t Þ; ∂t αðt Þ
0 o αðt Þ o 1; x A Ω; t A ð0; T Þ; ð21Þ
with zero initial condition uðx; 0Þ ¼ 0; x A Ω and under two different types of boundary conditions: (a) Full Dirichlet boundary conditions π x π x 1 2 uðx; t Þ ¼ t 2 sin sin ; x A Γ D ¼ ∂Ω; t A ð0; T Þm 2 2 ð22Þ (b) Mixed boundary conditions π x π x 1 2 sin ; uðx; t Þ ¼ t 2 sin 2 2 t A ð0; T Þ;
x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; ð23aÞ
∂uðx; t Þ ∂ t 2 sin ðπ x1 =2Þ sin ðπ x2 =2Þ ¼ ; ∂n ∂n
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ; ð23bÞ 2 αðt Þ 2 2 =Γ ð3 αðt ÞÞÞ þ ðπ t =2Þ sin ðπ x1 =2Þ where Q ðx; tÞ ¼ ð2t sin ðπ x2 =2Þ . The exact solution is
Fig. 13. Convergence rate (RMSE and RMSEx) of the present method with dt ¼ 0.01 by using different RBF formulations (11(b)) and (11(c)) in Example 3. (a) Full Dirichlet boundary conditions and (b) mixed boundary conditions.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
8
Table 1 Numerical errors (RMSE) by using the RBF formulation (11(b)) with both uniform and irregular nodes. Node distribution
RMSE
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
1.36E 04 4.50E 04
3.57E 03 6.38E 03
3.61E 04 1.90E 04
1.68E 03 2.53E 03
5.99E 03 6.17E 03
8.95E 03 1.07E 02
Table 2 Numerical errors (RMSEx) by using the RBF formulation (11(b)) with both uniform and irregular nodes. Node distribution
RMSEx
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
4.25E 03 4.72E 03
9.16E 03 3.15E 02
1.68E 02 1.38E 02
1.96E 02 6.23E 02
1.79E 02 1.61E 02
2.33E 02 6.18E 02
Table 3 Numerical errors (RMSE) by using the RBF formulation (11(c)) with both uniform and irregular nodes. Node distribution
RMSE
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
1.77E 04 1.89E 04
2.58E 03 3.69E 03
1.73E 04 1.62E 04
5.61E 04 7.94E 04
6.16E 03 6.17E 03
1.08E 02 1.10E 02
Table 4 Numerical errors (RMSEx) by using the RBF formulation (11(c)) with both uniform and irregular nodes. Node distribution
RMSEx
Uniform nodes Irregular nodes
Ex. 1(a)
Ex. 1(b)
Ex. 2(a)
Ex. 2(b)
Ex. 3(a)
Ex. 3(b)
2.47E 03 1.65E 03
6.73E 03 1.72E 02
7.28E 03 7.64E 03
8.67E 03 2.72E 02
1.31E 02 1.40E 02
1.61E 02 3.16E 02
different time steps (dt ¼0.02, 0.01) in Example 3. From Figs. 12 and 13, we can see that the present method with node density (dh ¼1/8, namely, M¼293) provides the acceptable numerical results (the error is less than 1%). And it has the similar conclusion with Examples 1 and 2 that the numerical accuracy of Example 3 (a) is better than those of Example 3(b) by using the present method with different basis functions (11(b))–(11(c)).
1 α =0.3 α =0.5 α =0.7 α =0.9
u(-0.5,0,t)
0.8 0.6
Next we place irregular nodes in the computational domain Ω1 as shown in Fig. 13 to study the meshless feature of RBF-type methods. Tables 1–4 list the present numerical errors (RMSE and RMSEx) with the RBF formulations (11(b))–(11(c)) by using both uniform and irregular nodes with node number M¼484 and time step dt¼0.01. From Tables 1–4, it can be seen that the present method works equally well with irregular nodes. And the present method with RBF formulations (11(c)) can perform better than RBF formulation (11 (b)), especially for the partial derivative term ð∂uðxi ; TÞÞ=∂x1 .
0.4 0.2 0
0
0.5
1 t
1.5
2
Fig. 14. Temporal evolution of uð 0:5; 0; t Þ with different time fractional derivative order α by using the present method with RBF formulations (11(c)).
uðx; t Þ ¼ t 2 sin
π x 1
2
sin
π x 2
2
;
x A Ω; t A ð0; T Þ:
ð24Þ
Uniformly distributed nodes in Fig. 1(a) are used. Figs. 12 and 13 show the convergence rate (RMSE and RMSEx) of the present method with different basis functions (11(b))–(11(c)) by using
Example 4. The last case is the following fractional diffusion model in Ω1 ∂α uðx; t Þ ¼ Δuðx; t Þ; ∂t α
0 o α o 1; x A Ω; t A ð0; T Þ;
ð25Þ
with initial condition uðx; 0Þ ¼ f ðxÞ; x A Ω and under the mixed boundary conditions
uðx; t Þ ¼ 0; x A Γ D ¼ ðx1 ; x2 Þjx1 ¼ 0; 2 ; t A ð0; T Þ; ð26aÞ
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
∂uðx; t Þ ¼ 0; ∂n where ( f ðxÞ ¼
x A Γ N ¼ ðx1 ; x2 Þjx2 ¼ 0; 2 ; t A ð0; T Þ;
2x1 þ2;
1 r x1 r1=2
ð2 2x1 Þ=3;
1=2 rx1 r 1
ð26bÞ
:
Finally, we set dt ¼0.01 and place uniform nodes (M¼488) in the computational domain Ω1 as shown in Fig. 1. Fig. 14 plots the temporal evolution of uð 0:5; 0; t Þ with different time fractional derivative order α by using the present method with RBF formulations (11(c)). Numerical investigation shows that the concentration u decays very fast at the beginning of diffusion process and it becomes more and more slowly at late times. This anomalous diffusion phenomenon becomes more and more obviously with the smaller time fractional derivative order α. 4. Conclusions This paper applies the method of approximate particular solutions (MAPS) to constant- and variable-order time fractional diffusion models. Discretization of the time-fractional and spatial derivatives is respectively done by a finite difference scheme and the MAPS. The present numerical experiments verify that the proposed method is a competitive meshless collocation method for constant- and variable-order time fractional diffusion models, and it provides satisfactory solutions and converges to the exact solutions with the increasing interpolation node number under
9
Moreover, it should be mentioned that the present numerical solutions of variable-order time fractional diffusion models are more sensitive to the time discretization formulation and depend on the function type of the fractional derivative order. Therefore, coupling MAPS with other higher order time discretization formulations [49,50] will be a sensible way to solve variable-order time fractional diffusion models. We leave the detailed numerical simulations to our future studies.
Acknowledgments The work described in this paper was supported by the National Science Funds of China (Grant nos. 11302069 and 11372097), the Fundamental Research Funds for the Central Universities, Hohai University (Grant no. 2013B32814), the National Basic Research Program of China (973 Project no. 2010CB832702), the National Science Funds for Distinguished Young Scholars of China (Grant no. 11125208), the 111 Project (Grant no. B12032), a CERG Grant of the Hong Kong Research Grant Council, and a FRG Grant of Hong Kong Baptist University.
Appendix A. Time discretization of variable-order time fractional derivative term The time discretization formulation (6) of variable-order time fractional derivative term (4) is derived as follows:
Z t ∂u x; η ∂αðt Þ uðx; t Þ 1 dη ¼ αðtÞ ∂η Γ ð1 α ðt ÞÞ 0 ∂t αðtÞ t η Z ðj þ 1Þτ k ∂u x; ξ 1 dξ ∑ α ðt j þ 1 Þ ∂ξ Γ 1 α t j þ 1 t j ¼ 0 jτ kþ1 ξ Z ðj þ 1Þτ k u x; t 1 dξ j þ 1 u x; t j ¼ ∑ τ Γ 1 α t j þ 1 t k þ 1 ξαðtj þ 1 Þ jτ j¼0 Z ðk j þ 1Þτ k u x; t 1 dη j þ 1 u x; t j ¼ ∑ τ Γ 1 α t k j þ 1 ηαðtk j þ 1 Þ ðk jÞτ j¼0 Z ðj þ 1Þτ k u x; t 1 dη k j þ 1 u x; t k j ¼ ∑ τ Γ 1 α t j þ 1 ηαðtj þ 1 Þ jτ j¼0 Z ðj þ 1Þτ k u x; t 1 dη k j þ 1 u x; t k j ¼ ∑ τ Γ 1 α t j þ 1 jτ ηαðtj þ 1 Þ j¼0
i h τ αðtj þ 1 Þ u x; t k j þ 1 u x; t k j ðj þ 1Þ1 αðt j þ 1 Þ j1 αðt j þ 1 Þ ; j ¼ 0 Γ 2 α tj þ 1 8 i α ðt j þ 1 Þ k > h > > τ αðt1 Þ uk þ 1 uk þ ∑ τ u x; t k j þ 1 u x; t k j ðj þ 1Þ1 αðtj þ 1 Þ j1 αðt j þ 1 Þ ; < Γ ð2 αðt 1 ÞÞ j ¼ 1Γ 2 α t j þ 1 ¼ > > > τ αðt1 Þ u1 u0 ; k ¼ 0 : Γ ð2 αðt 1 ÞÞ 8 k > > < a0 uk þ 1 uk þ ∑ aj bj uk j þ 1 uk j ; k Z 1 j¼1 ¼ > > : a0 u1 u0 ; k ¼ 0 k
¼ ∑
both uniform and irregular distributions. Numerical investigations show that the present method with RBF formulation (11(c)) performs better than that with RBF formulations (11(a)) and (11 (b)), and it cannot provide the correct numerical results with the partial derivative term ðð∂uðx; t Þ=∂nÞ or ð∂uðx; t Þ=∂x1 ÞÞ by using RBF formulation (11(a)) for 2D fractional diffusion equations.
kZ1
where a0 ¼
τ αðt1 Þ τ α ðt j þ 1 Þ ; ; a ¼ Γ ð2 α ðt 1 ÞÞ j Γ 2 α t j þ 1
h i bj ¼ ðj þ 1Þ1 αðtj þ 1 Þ j1 αðtj þ 1 Þ ;
j ¼ 1; 2; …; k
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i
10
Z.-J. Fu et al. / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
References [1] Hilfer R. Applications of fractional calculus in physics. Singapore: World Scientific; 2000. [2] Metzler R, Klafter J. The random walk's guide to anomalous diffusion: a fractional dynamics approach. Phys Rep 2000;339(1):1–77. [3] Wyss W. The fractional diffusion equation. J Math Phys 1986;27:2782–5. [4] Gorenflo R, Luchko Y, Mainardi F. Wright functions as scale-invariant solutions of the diffusion-wave equation. J Comput Appl Math 2000;118(1–2):175–91. [5] Yang Q, Turner I, Liu F, Ilis M. Novel numerical methods for solving the timespace fractional diffusion equation in two dimensions. SIAM J Sci Comput 2011;33(3):1159–80. [6] Gorenflo R, Mainardi F, Moretti D, Pagnini G, Paradisi P. Discrete random walk models for space-time fractional diffusion. Chem Phys 2007;284(1–2):521–41. [7] Fu Z-J, Chen W, Yang H-T. Boundary particle method for Laplace transformed time fractional diffusion equations. J Comput Phys 2013;235:52–66. [8] Chechkin AV, Gorenflo R, Sokolov IM. Fractional diffusion in inhomogeneous media. J Phys A: Math Gen 2005;38(42):L679. [9] Sun H, Chen W, Chen Y. Variable-order fractional differential operators in anomalous diffusion modeling. Phys A: Stat Mech Appl 2009;388(21):4586–92. [10] Lorenzo CF, Hartley TT. Variable order and distributed order fractional operators. Nonlinear Dyn 2002;29(1–4):57–98. [11] Sun H, Chen W, Wei H. A comparative study of constant-order and variableorder fractional models in characterizing memory property of systems. Eur Phys J Spec Top 2011;193(1):185–92. [12] Tadjeran C, Meerschaert MM. A second-order accurate numerical method for the two-dimensional fractional diffusion equation. J Comput Phys 2007;220(2):813–23. [13] Murio DA. Implicit finite difference approximation for time fractional diffusion equations. Comput Math Appl 2008;56(4):1138–45. [14] Wang H, Wang K, Sircar T. A direct O(Nlog2N) finite difference method for fractional diffusion equations. J Comput Phys 2010;229(21):8095–104. [15] Ren J, Sun Z, Zhao X. Compact difference scheme for the fractional subdiffusion equation with Neumann boundary conditions. J Comput Phys 2013;232(1):456–67. [16] Zhao X, Sun Z. A box-type scheme for fractional sub-diffusion equation with Neumann boundary conditions. J Comput Phys 2011;230(15):6061–74. [17] Sun H, Chen W, Li C, Chen Y. Finite difference schemes for variable-order time fractional diffusion equation. Int J Bifurc Chaos 2012;22(04):1250085. [18] Shen S, Liu F, Chen J, Turner I, Anh V. Numerical techniques for the variable order time fractional diffusion equation. Appl Math Comput 2012;218(22):10861–70. [19] Chen CM, Liu F, Anh V, Turner I. Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation. SIAM J Sci Comput 2010;32(4):1740–60. [20] Lin R, Liu F, Anh V, Turner I. Stability and convergence of a new explicit finitedifference approximation for the variable-order nonlinear fractional diffusion equation. Appl Math Comput 2009;212(2):435–45. [21] Langlands TAM, Henry BI. The accuracy and stability of an implicit solution method for the fractional diffusion equation. J Comput Phys 2005;205(2): 719–36. [22] Cui M. Compact finite difference method for the fractional diffusion equation. J Comput Phys 2009;228(20):7792–804. [23] Zhuang P, Liu F. Implicit difference approximation for the time fractional diffusion equation. J Appl Math Comput 2006;22(3):87–99. [24] Zhang Y, Sun Z, Zhao X. Compact alternating direction implicit scheme for the two-dimensional fractional diffusion-wave equation. SIAM J Numer Anal 2012;50(3):1535–55. [25] Chen C-M, Liu F, Turner I, Anh V. A Fourier method for the fractional diffusion equation describing sub-diffusion. J Comput Phys 2007;227(2):886–97. [26] Lin Y, Xu C. Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys 2007;225(2):1533–52.
[27] Jiang Y, Ma J. High-order finite element methods for time-fractional partial differential equations. J Comput Appl Math 2011;235(11):3285–90. [28] Li C, Zhao Z, Chen Y. Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion. Comput Math Appl 2011;62(3):855–75. [29] Fix GJ, Roof JP. Least squares finite-element solution of a fractional order twopoint boundary value problem. Comput Math Appl 2004;48(7–8):1017–33. [30] Katsikadelis JT. The BEM for numerical solution of partial fractional differential equations. Comput Math Appl 2011;62(3):891–901. [31] Brunner H, Ling L, Yamamoto M. Numerical simulations of 2D fractional subdiffusion problems. J Comput Phys 2010;229(18):6613–22. [32] Chen W, Ye L, Sun H. Fractional diffusion equations by the Kansa method. Comput Math Appl 2010;59(5):1614–20. [33] Liu Q, Gu Y, Zhuang P, Liu F, Nie Y. An implicit RBF meshless approach for time fractional diffusion equations. Comput Mech 2011;48(1):1–12. [34] Chen CS, Fan CM, Wen PH. The method of approximate particular solutions for solving elliptic problems with variable coefficients. Int J Comput Methods 2011;8(03):545–59. [35] Ling L, Chen CS, Kwok TO. Adaptive method of particular solution For solving 3d inhomogeneous elliptic equations. Int J Comput Methods 2010;07(03):499–511. [36] Chen W, Fu ZJ, Chen CS. Recent advances on radial basis function collocation methods. Berlin: Springer; 2013. [37] Tsai C-C, Chen CS, Hsu T-W. The method of particular solutions for solving axisymmetric polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2009;33(12):1396–402. [38] Chen CS, Fan CM, Wen PH. The method of approximate particular solutions for solving certain partial differential equations. Numer Methods Partial Diff Equ 2012;28(2):506–22. [39] Zhu H. The method of approximate particular solutions for solving anisotropic elliptic problems. Eng Anal Bound Elem 2014;40:123–7. [40] Reutskiy SY. Method of particular solutions for nonlinear Poisson-type equations in irregular domains. Eng Anal Bound Elem 2013;37(2):401–8. [41] Kuo LH, Gu MH, Young DL, Lin CY. Domain type kernel-based meshless methods for solving wave equations. CMC-Comput Mater Contin 2013;33(3): 213–28. [42] Bustamante CA, Power H, Florez WF, Hang CY. The global approximate particular solution meshless method for two-dimensional linear elasticity problems. Int J Comput Math 2013;90(5):978–93. [43] Bustamante CA, Power H, Sua YH, Florez WF. A global meshless collocation particular solution method (integrated radial basis function) for twodimensional Stokes flow problems. Appl Math Modell 2013;37(6):4538–47. [44] Jiang T, Li M, Chen CS. The method of particular solutions for solving inverse problems of a nonhomogeneous convection-diffusion equation with variable coefficients. Numer Heat Transf, Part A: Appl 2012;61(5):338–52. [45] Yao G, Siraj ul I, Sarler B. Assessment of global and local meshless methods based on collocation with radial basis functions for parabolic partial differential equations in three dimensions. Eng Anal Bound Elem 2012;36(11): 1640–8. [46] Reutskiy SY. Method of particular solutions for solving PDEs of the second and fourth orders with variable coefficients. Eng Anal Bound Elem 2013;37(10): 1305–10. [47] Bustamante CA, Power H, Florez WF. A global meshless collocation particular solution method for solving the two-dimensional Navier–Stokes system of equations. Comput Math Appl 2013;65(12):1939–55. [48] Podlubny I. Fractional differential equations. San Diego: Academic Press; 1999. [49] Cuesta E, Lubich C, Palencia C. Convolution quadrature time discretization of fractional diffusion-wave equations. Math Comput 2006;75(254):673–96. [50] Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes. AIAA paper; 1981-1259.
Please cite this article as: Fu Z-J, et al. Method of approximate particular solutions for constant- and variable-order fractional diffusion models. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.09.003i