Boundary particle method for Laplace transformed ... - Semantic Scholar
Journal of Computational Physics 235 (2013) 52–66
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Boundary particle method for Laplace transformed time fractional diffusion equations Zhuo-Jia Fu a,b, Wen Chen a,b,⇑, Hai-Tian Yang b a b
College of Mechanics and Materials and College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing, Jiangsu 210098, PR China State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China
a r t i c l e
i n f o
Article history: Received 19 March 2012 Received in revised form 1 October 2012 Accepted 18 October 2012 Available online 6 November 2012 Keywords: Boundary particle method Laplace transform Numerical inverse Laplace transform Meshless Time fractional derivative Anomalous diffusion
a b s t r a c t This paper develops a novel boundary meshless approach, Laplace transformed boundary particle method (LTBPM), for numerical modeling of time fractional diffusion equations. It implements Laplace transform technique to obtain the corresponding time-independent inhomogeneous equation in Laplace space and then employs a truly boundary-only meshless boundary particle method (BPM) to solve this Laplace-transformed problem. Unlike the other boundary discretization methods, the BPM does not require any inner nodes, since the recursive composite multiple reciprocity technique (RC-MRM) is used to convert the inhomogeneous problem into the higher-order homogeneous problem. Finally, the Stehfest numerical inverse Laplace transform (NILT) is implemented to retrieve the numerical solutions of time fractional diffusion equations from the corresponding BPM solutions. In comparison with finite difference discretization, the LTBPM introduces Laplace transform and Stehfest NILT algorithm to deal with time fractional derivative term, which evades costly convolution integral calculation in time fractional derivation approximation and avoids the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. Error analysis and numerical experiments demonstrate that the present LTBPM is highly accurate and computationally efficient for 2D and 3D time fractional diffusion equations. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Fractional diffusion equation [1] is considered a recent alternative model to describe anomalous diffusion phenomena [2] in a wide range of engineering and physics fields [3], such as electron transportation [4], seepage [5], magnetic plasma [6], dissipation [7] and turbulence [8], just to mention a few. Compared with the nonlinear models, the fractional anomalous diffusion models require fewer and easily available parameters to accurately describe the physical phenomena. Like the nonlinear models, the analytical solutions of fractional derivative models are largely not available, except for the cases with simple initial and boundary conditions by using Laplace transform [1,9] and Fourier–Laplace transforms. Hence the numerical simulation [10–12] plays an essential role in the analysis of fractional diffusion models. Nowadays the finite difference method (FDM) [13–15] is a popular and dominant numerical technique for time and spatial discretization of fractional derivative equations. Its convergence, accuracy, and stability have extensively been discussed in the literatures [16–19]. However, due to the memory (non-local) effect arising from the convolution integral expressions of time (space) fractional derivative term, the FDMs in these cases are no longer of a finite local discretization scheme and ⇑ Corresponding author at: College of Mechanics and Materials and College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing, Jiangsu 210098, PR China. Tel.: +86 25 83786873; fax: +86 25 83736860. E-mail address: [email protected] (W. Chen). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.10.018
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
53
encounter an exponential increase of computing costs with advancing time and large computational domain. Consequently, numerical results in most reports are restricted in 1D and 2D problems. This study considers time fractional derivative equations, which only has fractional derivative in time and integer differential operator in space. Most reports employ the finite difference method for temporal discretization but introduce the other numerical methods for spatial discretization, such as Fourier method [20], spectral method [21], finite element method [22–25], boundary element method [26], and radial basis function meshless collocation method [27–30]. In comparison with traditional FDMs, these methods can mitigate, to a certain extent, computing costs for large computational domain problems. But most of these methods still employ finite difference schemes for temporal discretization, and inevitably encounter an exponential increase of computing costs with advancing time and thus have low efficiency in the simulation of long timehistory fractional diffusion. In order to overcome the drawback of FDM for temporal discretization, this paper presents a novel approach, called Laplace transformed boundary particle method (LTBPM), for efficient modeling of time fractional diffusion equations. The LTBPM introduces Laplace transform to convert time fractional diffusion equation into a time-independent inhomogeneous equation in Laplace space. To retrieve the final time-dependent solutions from Laplace space solutions, the present LTBPM chooses a sophisticated Stehfest numerical inverse Laplace transform (NILT) [31], which is originated from Gaver [32] and applied to various scientific and engineering problems. It should be mentioned that, like the other NILT algorithms, the Stehfest algorithm is also an ill-posed problem and the magnification of truncation error results in a loss of numerical accuracy. For this reason, it is highly important to employ a highly-accurate and efficient spatial discretization method to solve the transformed time-independent diffusion equation in Laplace space. Fortunately, the authors recently develop boundary particle method (BPM) [33–35], which only requires boundary-only meshless discretization of inhomogeneous problem with integration-free merit and is very efficient and accurate for high-dimensional problems. In the present method, the BPM is used to solve the corresponding time-independent inhomogeneous equation in Laplace space. Therefore, the present LTBPM employs Laplace transform and Stehfest numerical inverse Laplace transform instead of the step-by-step FDMs for temporal discretization, which evades costly convolution integral calculation in time fractional derivative approximation and avoids the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. The rest of this paper is organized as follows. Section 2 introduces Laplace transformed boundary particle method (LTBPM) for numerical modeling of time fractional diffusion equations, followed by numerical results and discussions in Section 3. Some remarks are drawn upon the numerical results in Section 4.
2. Methodology 2.1. Time fractional diffusion problem Without loss of generality, we consider a time fractional diffusion equations in a bounded domain X with piecewise smooth boundary @ X ¼ CD þ CN ðCD \ CN ¼ ;Þ;
@ a uðx; tÞ ¼ DDuðx; tÞ þ ~ v ruðx; tÞ kuðx; tÞ þ Q ðx; tÞ; @t a
0 < a < 1;
x 2 X; t 2 ð0; TÞ;
ð1Þ
with boundary conditions
uðx; tÞ ¼ g 1 ðx; tÞ;
x 2 CD ;
@uðx; tÞ ¼ g 2 ðx; tÞ; @n
x 2 CN ;
t 2 ð0; TÞ;
ð2aÞ
t 2 ð0; TÞ;
ð2bÞ
and initial condition
uðx; 0Þ ¼ u0 ðxÞ;
x 2 X;
ð3Þ
where Q ðx; tÞ; g 1 ðx; tÞ; g 2 ðx; tÞ and u0 ðxÞ are known functions; x ¼ ðx; yÞ for 2D problem and x ¼ ðx; y; zÞ for 3D problem; D the diffusion coefficient, k the reaction coefficient, ~ v velocity vector, n the unit outward normal, T the total time to be considered, a and @t@ a the Caputo fractional derivative of order a with respect to t defined by
@ a uðx; tÞ 1 ¼ @t a Cð1 aÞ
Z 0
t
@uðx; gÞ dg ; @ g ðt gÞa
for more details see Podlubny [1].
0 < a < 1;
ð4Þ
54
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 1. The roadmaps of the LTBPM and the FDM for time fractional diffusion problems.
2.2. Implementation procedure Fig. 1 shows the roadmaps of the LTBPM and the FDM for time fractional diffusion problems. Unlike the one-step procedure in the FDM, the LTBPM includes three-step: (a) Laplace transform, (b) Boundary particle method (BPM), and (c) Numerical inverse Laplace transform (NILT). 2.2.1. Step (a) Laplace transform Laplace transform converts the time fractional diffusion problems from time domain to Laplace space domain. Let f(t) be a causal time domain function as an independent variable t P 0 and ~f ðsÞ denotes its image in Laplace space domain. Laplace transform is given by
~f ðsÞ ¼ Lðf ðtÞÞ ¼
Z
1
f ðtÞest dt;
ð5Þ
0
where s is Laplace transform parameter, and unless otherwise specified the quantities in Laplace space domain are denoted by an over tilde. Laplace transform of the Caputo fractional derivative [1] can be written as
a @ f ðtÞ ¼ sa ~f ðsÞ sa1 f ð0Þ; L a @t
0 < a < 1;
t P 0:
ð6Þ
Applying Laplace transform to Eqs. (1) and (2) produces
~ ðx; sÞ ¼ Fðx; sÞ; Ru ~ ðx; sÞ ¼ g~1 ðx; sÞ; u ~ ðx; sÞ @u ¼ g~2 ðx; sÞ; @n
0 < a < 1;
x 2 X;
x 2 CD ; x 2 CN ;
ð7Þ ð8aÞ ð8bÞ
where R ¼ DD þ ~ v r ðk þ sa Þ; Fðx; sÞ ¼ Q~ ðx; sÞ sa1 u0 ðxÞ. It is of interest to point out that, in comparison with finite difference discretization, Laplace transform evades costly convolution integral calculation of time fractional derivative and avoids the effect of time step on numerical accuracy and stability. 2.2.2. Step (b) Boundary particle method In the Laplace space, the boundary particle method (BPM) is implemented to solve Laplace-transformed time-indepen~ can be expressed as dent inhomogeneous problems (7) and (8) via meshless boundary nodes. The corresponding solution u
~p ; ~¼u ~h þ u u
ð9Þ
~ h and u ~ p are the homogeneous and the particular solutions, respectively. The particular solution u ~ p satisfies where u
~ p g ¼ F; Rfu
ð10Þ
~ h has to satisfy not only the but does not necessarily satisfy boundary conditions. In contrast, the homogeneous solution u corresponding homogeneous equation
~ h g ¼ 0; Rfu
ð11Þ
but also boundary conditions
~p ; ~ h jC ¼ g~1 u u D
ð12aÞ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~h ~ @u ¼ g~2 @ up : @n CN @n
55
ð12bÞ
~ h of Eqs. (11) and (12) can be efficiently calculated by boundary-type In the BPM procedure, The homogeneous solution u numerical techniques with the nonsingular general solutions [36], Trefftz functions [37] and singular fundamental solutions [38]. ~ p , the BPM employs the recursive composite multiple reciprocity method (RC-MRM) To evaluate the particular solution u [33] to annihilate the frequently encountered source term Fðx; sÞ in Eq. (10) by using composite differential operators
Lm . . . L2 L1 fFðx; sÞg ffi 0;
ð13Þ
where L1, L2, . . . Lm are differential operators of the same or different types, and Fðx; sÞ represents the polynomial, exponential and trigonometric functions or the combination of these functions in this study, which frequently appear in engineering problems. Consequently, Eq. (10) is equivalently transformed into the following high-order homogeneous problem
8 ~ p ðx; sÞ ¼ 0 x 2 X Lm . . . L2 L1 Ru > > > > .. < . ; > > L1 Ru ~ p ðx; sÞ ¼ L1 ðFðx; sÞÞ x 2 @ X > > : ~ p ðx; sÞ ¼ Fðx; sÞ Ru x 2 @X
0 < a < 1:
ð14Þ
Rearrange Eqs. (11) (12) and (14), we obtain the following high-order homogeneous problem
8 ~ ðx; sÞ ¼ 0 Lm . . . L2 L1 Ru > > > > > .. > > . > > > < ~ L1 Ruðx; sÞ ¼ L1 ðFðx; sÞÞ > ~ ðx; sÞ ¼ Fðx; sÞ > Ru > > > > ~ ðx; sÞ ¼ g~1 ðx; sÞ; > u > > > ~ ðx;sÞ : @u ¼ g~2 ðx; sÞ; @n
x2X
x 2 @X x 2 @X
0 < a < 1:
ð15Þ
x 2 CD x 2 CN
Then the solution of the transformed inhomogeneous Eqs. (7) and (8) is equal to that of the homogeneous problem (15). can be expressed as a linear combination of singular fundamental solutions with unHence, the approximate solution u known coefficients faij g as follows
ðxÞ ¼ u
m X N X aij Gi ðx yj Þ;
x 2 X;
ð16Þ
i¼0 j¼1
where N represents the number of source points, fyj g denotes the source points placed on the fictitious boundary outside the physical domain X, and Gi means the high-order singular fundamental solution of the differential operator Li (set L0 ¼ R). Some frequently-used high-order singular fundamental solutions are listed in Appendix. It is of interest to point out that the BPM is a truly boundary-only meshless and integration-free numerical technique. 2.2.3. Step (c) Numerical inverse Laplace transform Numerical inverse Laplace transform (NILT) is implemented to convert the BPM solutions from Laplace space domain to time domain. As the NILT is an ill-posed problem, the truncation error magnification appears and results in a loss of numerical accuracy. Here we employs the well-established Stehfest NILT algorithm originated from Gaver [32] to retrieve the timedependent solution of the original problem (1)–(3) from the BPM solutions in Laplace space. An approximate value fs of the Stehfest NILT of f(t) for a specific time t is given by
fs ðtÞ ¼
M ln 2 X ln 2 i ; V i ~f N t i¼1 t
ð17Þ
where ~f N is the Laplace space solutions, and M þi 2
V i ¼ ð1Þ
min ði;M 2Þ
X
k¼½iþ1 2
M
M 2
k 2 ð2kÞ! ; k !k!ðk 1Þ!ði kÞ!ð2k iÞ!
in which [C] denotes the nearest integers less than or equal to C. For a specific instant t, one should solve M boundary value problems (7) and (8) for the corresponding Laplace transform parameter s ¼ lnt 2 i; i ¼ 1; 2; . . . ; M. Then the solution of the time fractional diffusion Eqs. (1)–(3) can be retrieved by using the expression (17).
56
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
2.3. Error analysis This section presents error analysis of the present LTBPM. Since all of Laplace transforms are analytical in Step (a), and LTBPM’s numerical error arises largely in Step (b) and Step (c). First, we discuss the numerical error in Step (c), i.e. Stehfest numerical inverse Laplace transform (NILT). Abate et al. [39,40] did extensive numerical experiments to investigate the parameter effect on numerical accuracy, and summarized the following conclusion: ‘‘If p significant digits are desired, then let M be the positive integer [2.2p]; Given M, set the system precision at q ¼ ½1:1M; Given the transform ~f N and the argument t, then calculate fs by Eq. (17)’’. According to Abate’s conclusion, we obtain the following heuristic error estimation: ~
~
Remark 1. If the error of input data is 10ðqþ1Þ 6 kf N~f k 6 10q with even positive integer M via q ¼ ½1:1M, then the final kf k s f k error is 10ðpþ1Þ 6 kfkf 6 10p , where M ¼ ½2:2p. k It can be found from Remark 1 that the Stehfest algorithm tends to demand highly precise input data ~f N in order to yield good accuracy in the following NILT calculations. For example, if input data ~f N is single-precision data, namely, q = 7 with M = 8, then p = 4. If input data ~f N are double-precision data, namely, q = 16 with M = 16, then p = 7.In the present methodology, the input data ~f are obtained by solving time-independent inhomogeneous problems (7) and (8) in Laplace space. Therefore, it is necessary to employ a highly accurate numerical method in Step (b) to guarantee high accuracy in the final solutions. For this purpose, we implement the boundary particle method (BPM) with singular fundamental solutions in Step (b) to obtain highly accurate solution in Laplace space. As we will see in the following, the Laplace-transformed governing Eq. (7) degenerates to modified Helmholtz equation when the parameter ~ v ¼ ½0; 0. In BPM’s error analysis, we assume that the solution u~, the source term function F, and the ~1 are smooth enough on X. At first, the following homogeneous modified Helmholtz equation boundary condition function g with Dirichlet boundary condition under a unit disc is considered
~ ¼ 0; ðD k2 Þu
ð18Þ
~ jC ¼ g~1 : u D
ð19Þ
According to Theorem 2.3 in [41], we obtain a simplified error estimation N
~u kL1 ðXÞ 6 skg~1 kL1 ðC Þ ku D
d
N
1d
ð20Þ
;
is the corresponding approximation solution, the fictitious boundary parameter d > 1 represents the radius of a conwhere u centric disc Xd with unit disc X, N denotes the number of source points on @ Xd or collocation points on CD, s is a constant independent of N. The next step turns to consider the following inhomogeneous modified Helmholtz equation
~ ¼ F; ðD k2 Þu
ð21Þ
~ jC ¼ g~1 ; u D
ð22Þ
in which the source term F can be eliminated by using second-order differential operator L1, namely, L1 fFg ¼ 0 on X. ~ p and u ¼u h þ u p be the exact and numerical solution of the inhomogeneous problem (21) and (22), ~¼u ~h þ u Theorem 1. Let u Then p is the particular solution of Eq. (21), and ðD k2 Þu p ¼ F. where u N
~u kL1 ðXÞ 6 s1 kg~1 u p k L1 ðC Þ ku D where
d
N
1d
þ eF ;
ð23Þ
s1 is a constant independent of N, and eF satisfies
eF
8 > s2 kFkL1 ðCD Þ dN ; L1 ¼ D; > > < N ¼ s3 kFkL1 ðCD Þ d N ; L1 ¼ D k21 ; ; 1d > > > : s kFk 1 dN ; L1 ¼ D þ k21 4 L ðCD Þ
ð24Þ
in which s2, s3 and s4 are the constant independent of N. ^ h be the exact solution of the following problem Proof. Let u
^ h ¼ 0; ðD k2 Þu
ð25Þ
p : ^ h jC ¼ g~1 u u D
ð26Þ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
57
^¼u ^h þ u p ; then we have Let u
~u kL1 ðXÞ 6 ku ~u ^ kL1 ðXÞ þ ku ^u kL1 ðXÞ : ku
ð27Þ
By using Eq. (20), we have N
^u kL1 ðXÞ ¼ ku ^h u h kL1 ðXÞ 6 s1 kg~1 u p kL1 ðC Þ ku D
d
N
1d
ð28Þ
:
Notice that in X
~u ^Þ ¼ F F; ðD k2 Þðu
ð29Þ
and on CD
~u ^ÞjC ¼ g~1 ðg~1 u p þ u p Þ ¼ 0: ðu D
ð30Þ
According to a priori estimate [42], we obtain
1 : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ
ð31Þ
Consider the following problems
L1 v ¼ 0;
ð32Þ
v jC
ð33Þ
D
¼ F;
where F and F are its exact and numerical solution, respectively. For L1 ¼ D k21 , Eq. (31) can be rewritten by using Eq. (20)
1 6 s3 kFk 1 ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ
d
N N
1d
ð34aÞ
:
For L1 ¼ D, Eq. (31) can be rewritten by using Theorem 2 in [43]
1 6 s2 kFk 1 dN : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ For L1 ¼ D þ
k21 ,
ð34bÞ 2
Eq. (31) can be rewritten by using Theorem 3 with q > R in [44]
1 6 s4 kFk 1 dN : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ
ð34cÞ
Therefore, the conclusion of the theorem follows from (27), (28) and (34). It can be observed from Eq. (23) that the BPM with singular fundamental solution has nearly exponential convergence for some inhomogeneous modified Helmholtz problems. Hence, the BPM error decreases exponentially as N ! 1 and d ! 1. As a price pays for such a high accuracy, the condition number of BPM interpolation matrix increases exponentially with N ! 1 and d ! 1. Fix on d, the BPM error tends to decrease exponentially and then oscillate slightly below a certain error level with an increase of N, due to the machine round-off error which becomes dominant in the final error at large value of N. Extensive numerical experiments [33,34,37] also demonstrate this observation. h 3. Numerical results and discussion In this section, the efficiency, accuracy and convergence of the proposed LTBPM are tested to four time fractional diffusion equations. The LTBPM solutions are compared with domain-type radial basis function (DRBF) solution [29] and analytical solution. For quantitative examinations, the notations of the maximum absolute and relative average error norms are defined below
ðjÞj; MerrðuÞ ¼ max juðjÞ u 16j6NT
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPNT u j¼1 ðuðjÞ u ðjÞÞ2 ; RerrðuÞ ¼ t PNT 2 j¼1 ðuðjÞÞ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPNT u j¼1 ðu;x ðjÞ u ;x ðjÞÞ2 ; Rerrðu;x Þ ¼ t PNT 2 j¼1 ðu;x ðjÞÞ
ð35aÞ
ð35bÞ
ð35cÞ
ðjÞ; u ;x ðjÞ and uðjÞ; u;x ðjÞ are the analytical and numerical solutions at xj , respectively, and NT denotes the total number where u of uniform test points in the domain of interest. In this study, NT is taken to be 121 in 2D cases and 1331 in 3D cases. In order to investigate the BPM accuracy in Step (b), we define
58
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~ Þ ¼ max ðMerrðu ~ ÞÞ; merrðu
ð36aÞ
~ ÞÞ; ~ Þ ¼ max ðRerrðu rerrðu
ð36bÞ
~;x ÞÞ; ~ ;x Þ ¼ max ðRerrðu rerrðu
ð36cÞ
16j6M
16j6M
16j6M
~; u ~ ;x are the solutions of Eqs. (7) and (8), i.e. the NILT input data ~f N in Eq. (17). where u Unless otherwise specified, the fictitious boundary parameter d = 6. For 2D cases, the source points are placed on a fictitious circle centering at geometric origin op and the fictitious boundary parameter d is its radius. For 3D cases, the fictitious boundary parameter d is defined as
d¼
yj xj : xj op
ð37Þ
Example 1. We first consider the time fractional advection-diffusion equation with the Dirichlet boundary conditions in a unit square domain X ¼ ½0; 1 ½0; 1. In order to compare with the reference results [29], the parameters in Eqs. (1)–(3) are h 2a i 2 xþy chosen as D ¼ 1; ~ v ¼ ½0; 0; k ¼ 0, Q ðx; tÞ ¼ C2tð3 ; g 1 ðx; tÞ ¼ t2 exþy ; u0 ðxÞ ¼ 0; where a ¼ 0:85; op = (0.5, 0.5) and aÞ 2t e the analytical solution uðx; tÞ ¼ t2 exþy , in which x ¼ ðx; yÞ 2 X. 2 exþy ; By applying the Laplace transform, we have the operators and functions R ¼ D sa ; Fðx; sÞ ¼ s43 þ s3 a g~1 ðx; sÞ ¼ s23 exþy in transformed problems (7) and (8). It is noted that the source term Fðx; sÞ can be annihilated by employing the differential operator L ¼ D 2 in the BPM.Table 1 displays the numerical accuracy of input data ~f obtained by the BPM N
1
with different boundary-only discretizations at T = 1. Numerical results verify that the BPM errors tend to decrease rapidly and then oscillate slightly below a certain error level with node refinements, which is in accordance with the BPM error analysis in Section 2.3. Fig. 2 shows the LTBPM error variations with respect to M in Eq. (17) by using different boundary node spacings (Dh ¼ 0:2; 0:1; 0:05) in comparison with the final errors obtained by the NILT with exact input data ~f . From Fig. 2, it can be found that all of the final errors tend to decrease rapidly and then increase with an increase of M, and all the turning points of error curves depend on the accuracy of NILT input data. It should be mentioned that in this case the NILT uses the BPM solutions ~f N (Dh ¼ 0:2; 0:1; 0:05) as input data to obtain the similar accuracy with exact input data ~f when M 6 14. It can be observed from Table 1 and Fig. 2 that when M 6 14 the loss of numerical accuracy mainly occurs in the NILT process rather than in the BPM computation. It is because the BPM solutions with Dh ¼ 0:2 are accurate enough to be the NILT input data, and the ill-posed NILT process results in the loss of numerical accuracy. This phenomenon can be also found in the method of particular solution with NILT for time-dependent problems in reference [45]. Furthermore, Table 2 presents the LTBPM errors with M ¼ 10 at T = 1 in comparison with the DRBF method. Numerical results show that the proposed LTBPM provides very accurate results by boundary-only discretization with node spacing Dh ¼ 0:2. In contrast, according to the convergence rate formula [29], the DRBF method requires whole domain discretization with smaller node spacing Dh ¼ 0:1 and very small time step Dt ¼ 0:004 to obtain the similar accuracy as the LTBPM. Namely, to achieve the similar accuracy, the LTBPM requires solving 20 20 linear systems with ðm þ 1Þ M ¼ 20 times, while the DRBF method requires solving 100 100 linear systems with T=Dt ¼ 250 times. Roughly speaking, the present study shows that the DRBF method requires more computational resources than the LTBPM. Moreover, in the time fractional derivative approximation, the FDM encounters an exponential increase of computing costs with advancing time because its memory requirements can grow in a dramatic way. In contrast, the LTBPM uses Laplace transform to evade this enormous effect on memory demand. These sharp comparisons indicate that the LTBPM is far more efficient for a long time-history fractional diffusion simulation. Example 2. The next case is a three-dimensional time fractional diffusion equation with the Dirichlet boundary conditions 2a
in a cube domain X ¼ ½0; 2 ½0; 2 ½0; 2. The parameters are Q ðx; tÞ ¼ ½C2tð3aÞ þ t 2 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ þ 6t2 , D ¼ 1; ~ v ¼ ½0; 0; k ¼ 1, u0 ðxÞ ¼ 0, g 1 ðx; tÞ ¼ t ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ, where a ¼ 0:7, op = (1, 1, 1) and the analytical solution is uðx; tÞ ¼ t2 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ, in which x ¼ ðx; y; zÞ 2 X. As in the first case, applying the Laplace transform yields the operators and source functions R ¼ D ðk þ sa Þ, 2 12 ~ 2 Fðx; sÞ ¼ ðs23 þ s3 1 ðx; sÞ ¼ s3 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ in the transformed equaa Þ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ þ s3 , g tions (7) and (8). The source term Fðx; sÞ can be annihilated by employing the differential operators L2 ¼ L1 ¼ D. 2
Table 1 The BPM errors with different boundary-only discretizations at T = 1. BPM
~Þ merrðu
~Þ rerrðu
~ ;x Þ rerrðu
Dh = 0.2 Dh = 0.1 Dh = 0.05
4.115e9 6.573e13 7.603e13
4.114e11 7.768e15 2.186e14
5.479e10 4.113e14 1.954e14
59
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 2. The final errors (a) Rerr(u) and (b) Rerr(u,x) with respect to M by using NILT with different input data (BPM solutions ~f N with Dh ¼ 0:2; 0:1; 0:05 and analytical solutions ~f ).
Table 2 The errors obtained by LTBPM with boundary-only discretization (Dh = 0.2) at T = 1 in comparison with the DRBF method with domain discretization (Dh = 0.1).
LTBPM DRBF(Dt = 0.1) [29] DRBF(Dt = 0.004)a a
Merr(u)
Rerr(u)
Rerr(u,x)
Rerr(u,y)
4.135e4 1.234e2 3.046e4
5.596e5 2.073e3 5.116e5
5.596e5 6.864e3 1.694e4
5.596e5 6.864e3 1.694e4
The DRBF results with Dt = 0.004 were obtained from the convergence rate formula OðDt2a Þ in reference [29].
Table 3 lists the BPM results at various time instants T by different boundary-only discretizations on the surface of 3D cube domain. It can be seen from Table 3 that the BPM with the same boundary discretization achieves the similar accuracy on various time instants T, even with small or large time instants T (T = 0.01, T = 100). It indicates that the present LTBPM avoids the effect of time step on numerical accuracy and stability encountered in the finite difference methods, and effectively simulates the long time-history fractional diffusion systems. Fig. 3 presents the LTBPM error variations with respect to M in Eq. (17) by using different boundary node spacings (Dh ¼ 0:4; 0:2) in comparison with the resulting errors by using exact input data in the NILT, and the corresponding LTBPM results at z = 1 cross section. This case verifies that the proposed LTBPM works equally well and efficiently for 3D fractional diffusion equations. It can be also observed from Fig. 3 that all of the final errors tend to decrease rapidly and then increase with increasing M, and all the turning points of error curves depend on the accuracy of NILT input data. Furthermore, the NILT uses the BPM solutions ~f N (Dh ¼ 0:2; 0:1; 0:05) as input data to obtain the similar accuracy with exact input data ~f when M 6 14. It also indicates that the final error is generated mainly by the Stehfest numerical inverse Laplace transform rather than by the BPM computation. It should be mentioned that Remark 1 is the heuristic error estimation. Hence we can estimate the optimal choice M is around 16 by using double-precision data type in the present numerical experiments. Moreover, the BPM numerical error in Step (b) (BPM NILT input data) is larger than the machine round-off error (exact NILT input data). Therefore, the optimal choice M may be less than 16 according to Remark 1. In addition, the numerical results in Examples 1 and 2 also demonstrate that the above discussions of the optimal choice M are basically correct, where M = 14 appears the best choice for simulation. So we choose M = 14 in the following numerical computation.
Table 3 The BPM errors with different boundary-only discretizations at various time instants T. BPM
Dh = 0.4
Dh = 0.2
T = 0.01 T=1 T = 100 T = 0.01 T=1 T = 100
~Þ rerrðu
~ ;x Þ rerrðu
1.045e9 1.045e9 1.044e9 4.475e12 4.455e12 2.967e12
7.464e9 7.464e9 7.457e9 1.986e11 2.666e11 3.147e11
60
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 3. LTBPM error variations with respect to M on different time instants (a) T = 0.01, (c) T = 1, (e) T = 100; and LTBPM solution at z = 1 cross section on different time instants (b) T = 0.01, (d) T = 1, (f) T = 100.
Example 3. This case tests the proposed LTBPM approach to the time fractional diffusion equation with zero Dirichlet/mixed boundary conditions in a square domain X ¼ ½1; 1 ½1; 1 [27]. The parameters in Eqs. (1)–(3) are chosen as D ¼ 1; ~ v ¼ ½0; 0; k ¼ 0, Q ðx; tÞ ¼ 0; g 1 ðx; tÞ ¼ g 2 ðx; tÞ ¼ 0, op = (0, 0), in which x ¼ ðx; yÞ 2 X. For Dirichlet boundary conditions, the analytical solution associated with initial condition
u0 ðxÞ ¼ cos is given by
p p x cos y ; 2 2
ð38Þ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
61
Fig. 4. The maximum norm of the analytical solution with Dirichlet and mixed boundary conditions, respectively, for a ¼ 0:1; 0:2; . . . ; 0:9 within time period T 2 ½0; 1.
Fig. 5. The final errors obtained by the NILT with exact input data for Dirichlet and mixed boundary conditions, respectively, at a ¼ 0:1; 0:2; . . . ; 0:9 within time period T 2 ½0; 1.
p p 1 uðx; tÞ ¼ Ea p2 t a cos x cos y : 2 2 2
ð39Þ
For mixed boundary conditions having Dirichlet boundary on x ¼ 1 and x ¼ 1 and Neumann boundary on y ¼ 1 and y ¼ 1. The initial condition is
u0 ðxÞ ¼ cos
p p x sin y : 2 2
ð40Þ
Then we have the analytical solution
p p 1 uðx; tÞ ¼ Ea p2 t a cos x sin y : 2 2 2
ð41Þ
where the one-parameter Mittag–Leffler function is defined as
Ea ðzÞ ¼
1 X
zj
j¼0
Cðaj þ 1Þ
; a > 0:
ð42Þ
62
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~ Þ) with Fig. 6. LTBPM error profiles (Rerr(u)) with boundary-only discretizations (a) Dh = 0.2, (c) Dh = 0.1; and the corresponding BPM error profiles (rerrðu boundary-only discretizations (b) Dh = 0.2, (d) Dh = 0.1 for Dirichlet problem in Example 3.
As in the previous cases, we employ the Laplace transform to get the expression of R ¼ D sa ; Fðx; sÞ ¼ sa1 u0 ðxÞ; ~2 ðx; sÞ ¼ 0 in transformed Eqs. (7) and (8). The source term Fðx; sÞ can be annihilated by employing the differential g~1 ðx; sÞ ¼ g 2 operator L1 ¼ D þ p2 in the BPM. It should be stressed that the BPM must add a constraint condition to guarantee the uniqueness of the solution of the transformed high-order homogeneous equations. The simplest way is to add the approximate formulation (16) to satisfy the governing Eq. (7) in few additional inner points. In this case, only one additional inner point (0, 0) is enough to get accurate and stable solutions. Fig. 4 displays the maximum norm of the analytical solutions with Dirichlet and mixed boundary condition problems, respectively, for a ¼ 0:1; 0:2; ; 0:9 and time period T 2 ½0; 1. It can be found that the maximum norms drop very quickly at the beginning instants, especially for small a. As time advances, the solution varies far less after a short period. Fig. 5 shows the final errors obtained by the NILT with exact input data for Dirichlet and mixed boundary condition problems, respectively, at a ¼ 0:1; 0:2; ; 0:9 and time period T 2 ½0; 1. We can see from Fig. 5 that the errors decrease with an decrease of fractional order a. ~ Þ) by using boundary-only Fig. 6 shows the LTBPM error profiles (RerrðuÞ) and the corresponding BPM error profiles (rerrðu discretization with node spacing Dh ¼ 0:2 and Dh ¼ 0:1 for Dirichlet boundary condition case under a ¼ 0:1; 0:2; ; 0:9 within T 2 ½0; 1. From Fig. 6, one can observe that the LTBPM achieves acceptable accuracy merely with Dh ¼ 0:2 boundaryonly discretization. With a further node refinement (Dh ¼ 0:1), the LTBPM almost obtains the same accuracy as the NILT with exact input data ~f as shown in Fig. 5. This means that the BPM solutions ~f N with Dh ¼ 0:1 are accurate enough to be the NILT input data, and the loss of numerical accuracy is mainly generated in the ill-posed NILT process due to its error magnification.
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
63
~ Þ) with Fig. 7. LTBPM error profiles (Rerr(u)) with boundary-only discretizations (a) Dh = 0.2, (c) Dh = 0.1; and the corresponding BPM error profiles (rerrðu boundary-only discretizations (b) Dh = 0.2, (d) Dh = 0.1 for mixed boundary problem in Example 3.
~ Þ) by using Fig. 7 further illustrates the LTBPM error profiles (RerrðuÞ) and the corresponding BPM error profiles (rerrðu boundary-only discretization with node spacing Dh ¼ 0:2 and Dh ¼ 0:1 for mixed boundary condition case under a ¼ 0:1; 0:2; . . . ; 0:9 within T 2 ½0; 1. The similar observations can be made as in the previous Dirichlet problem. Example 4. The last case is concerned with the fractional subdiffusion–convection problem with mixed boundary conditions v ¼ ½0:005; 0; k ¼ 0, Q ðx; tÞ ¼ 0, in a square domain X ¼ ½1; 1 ½1; 1. The corresponding parameters are D ¼ 1; ~ 5
u0 ðxÞ ¼ ðxþ1Þ , op = (0, 0), The Dirichlet boundary conditions g 1 ðx; tÞ ¼ x þ 1 on x ¼ 1 and x ¼ 1 and the Neumann boundary 16 conditions g 2 ðx; tÞ ¼ 0 on y ¼ 1 and y ¼ 1, in which x ¼ ðx; yÞ 2 X. @ We can derive the operators and source functions of the Laplace transformed Eqs. (7) and (8), R ¼ D þ 0:005 @x sa ; ~ Fðx; sÞ ¼ sa1 u0 ðxÞ; g~1 ðx; sÞ ¼ xþ1 ; g ðx; sÞ ¼ 0. It is easily to find that the source term Fðx; sÞ can be annihilated by 2 s employing the differential operators L3 ¼ L2 ¼ L1 ¼ D in the BPM. Due to the symmetry of this problem, Fig. 8 shows LTBPM solutions with boundary-only discretization (Dh ¼ 0:2) at the cross section (parallel to the x1-axis) with varied fractional order a for every 1/32 stepping up from the lower-right towards the diagonal. The proposed LTBPM utilizes d ¼ 2:5 and coarse discretization (Dh ¼ 0:2) to produce the similar solution profiles as given in the reference method [27]. It can be observed from Fig. 8 that the effect of convection is significant in the case of a ¼ 0:9, and the inflection points can more clearly be seen at the initial instants. As the fractional order a decreases,
64
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 8. LTBPM solution (cross section) of the fractional subdiffusion-convection problem with varied order a of fractional derivative by using boundary-only discretization (Dh = 0.2).
the effect of convection is weakened. It should be noting that the numerical solutions for different a vary tremendously near the left endpoint x ¼ 1. The above numerical results favourably suggest that sensors should be placed somewhere in 1 < x < 0 instead of in 0 < x < 1 in the experiment design to capture the dramatic variance of the solution.
4. Conclusions This paper develops a Laplace transformed boundary particle method for numerical modeling of time fractional diffusion equations. In contrast to the standard FDMs, the proposed scheme implements Laplace transform and Stehfest numerical Laplace inversion to handle the time fractional derivative. It first uses Laplace transform to convert time fractional diffusion equation into a time-independent inhomogeneous equation in Laplace space. And then it uses a truly boundary-only meshless BPM to discretize the spatial derivatives of this Laplace-transformed inhomogeneous equation only with boundary discretization nodes. Finally, it chooses Stehfest numerical Laplace inversion to retrieve the numerical solutions of time fractional diffusion equations from the corresponding BPM solutions in Laplace space. In comparison with the finite difference discretization, the present LTBPM introduces Laplace transform technique to evade costly convolution integral calculation in time fractional derivation approximation and to avoid the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. On the other hand, the BPM makes the present strategy very attractive to high-dimensional problems. Error analysis and numerical experiments demonstrate that the LTBPM is a competitive numerical method in the solution of 2D and 3D time fractional diffusion equations with the merits of mathematically simple, easy-to-program, inherently meshless, highly accurate, boundary-only discretization and integration-free.
65
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Moreover, it is observed from numerical results that the Stehfest numerical inverse Laplace transform has the error magnification in ill-posed process results, which results in loss of numerical accuracy. To solve this ill-posed issue, the present LTBPM implements the BPM to produce the highly accurate NILT input data and guarantee that the NILT still achieves high accuracy in the final results. Furthermore, it is of interest to point out that the multiple precision arithmetic [39] becomes a popular and active research area in science and engineering community. It can be introduced to further overcome the loss of numerical accuracy in the ill-posed NILT process, which is now under intense study. The further results will be reported in a subsequent paper. Acknowledgements The authors thank Prof. CS Chen and the anonymous reviewers of this article for their very helpful comments and suggestions to significantly improve the academic quality of this article. This work was supported by National Basic Research Program of China (973 Project No. 2010CB832702), R&D Special Fund for Public Welfare Industry (Hydrodynamics, Grant No. 201101014), National Science Funds for Distinguished Young Scholars (Grant No. 11125208) and Programme of Introducing Talents of Discipline to Universities (111 project, Grant No. B12032). The first author would like to thank Hohai University Training Program for Excellent Doctoral Dissertation (Grant No. 2011B14814), Jiangsu Province Graduate Students Research and Innovation Plan (Grant No. CX10B_203Z) for financial support. Appendix A In the appendix, D denotes Laplacian operator, r the gradient operator, D the diffusivity coefficient, k a real number qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 j~ vj 2 l ¼ ð2D Þ þ kD , j foundation stiffness, and r
v and ~r the velocity vector and distance vector, known as the wave number, ~
the Euclidean norm between the point x and the origin, i.e. r ¼ kxk. Ym and Km represent the m-order Bessel and modified Bessel functions of the second kind, respectively. Note that the fundamental solution to a differential operator is not unique. For the Laplacian, a constant may be included. However, it is often omitted for it does not have a significant effect on the accuracy of the numerical results. As discussed in [38], the high-order fundamental solution of the Laplacian operator
Dmþ1 can be derived as follow ( r2m ðC m ln r Bm Þ; 2p /m ðxÞ ¼ F 1 r 2m1 ; ð2mÞ! 4p
x 2 R2 x 2 R3
ðA1Þ
;
Cm 1 where C 0 ¼ 1; B0 ¼ 0; C mþ1 ¼ 4ðmþ1Þ 2 ; Bmþ1 ¼ 4ðmþ1Þ2
Cm mþ1
þ Bm .
The high-order fundamental solution of Helmholtz-type operator ðD þ k2 Þmþ1 [33] is given by mþ1n=2 /m Y m1þn=2 ðkrÞ; F ðxÞ ¼ Am ðkrÞ
x 2 Rn :
ðA2Þ 2 mþ1
The high-order fundamental solution of modified Helmholtz-type operator ðD k Þ
/m F ðxÞ ¼ Am ðkrÞ
mþ1n=2
K m1þn=2 ðkrÞ;
[33] is given by
x 2 Rn :
ðA3Þ 2 mþ1
The high-order fundamental solution of modified convection-diffusion-type operator ðDD þ v r k Þ
[33] is given
by ~ v~r
mþ1n=2 /m e 2D K m1þn=2 ðlrÞ; x 2 Rn ; F ðxÞ ¼ Am ðlrÞ
ðA4Þ
2
where Am ¼ Am1 =2mk ; A0 ¼ 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Reports 339 (2000) 1–77. R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000. H. Scher, E.W. Montroll, Anomalous transit-time dispersion in amorphous solids, Physical Review B 12 (1975) 2455–2477. R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini, P. Paradisi, Discrete random walk models for space-time fractional diffusion, Chemical Physics 284 (2007) 521–541. D. del-Castillo-Negrete, B.A. Carreras, V.E. Lynch, Front dynamics in reaction-diffusion systems with levy flights: a fractional diffusion approach, Physical Review Letters 91 (2003) 018302. T.L. Szabo, J. Wu, A model for longitudinal and shear wave propagation in viscoelastic media, Journal of Acoustical Society of America 107 (2000) 2437– 2446. I.M. Sokolov, J. Klafter, A. Blumen, Ballistic versus diffusive pair dispersion in the Richardson regime, Physical Review E 61 (2000) 2717–2722. R. Gorenflo, Y. Luchko, F. Mainardi, Wright functions as scale-invariant solutions of the diffusion-wave equation, Journal of Computational and Applied Mathematics 118 (2000) 175–191. K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Computer Methods in Applied Mechanics and Engineering 194 (2005) 743–773. S. Momani, An algorithm for solving the fractional convection–diffusion equation with nonlinear source term, Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1283–1290.
66
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
[12] P.P. Valko, J. Abate, Numerical inversion of 2-D Laplace transforms applied to fractional diffusion equations, Applied Numerical Mathematics 53 (2005) 73–88. [13] C. Tadjeran, M.M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, Journal of Computational Physics 220 (2007) 813–823. [14] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equations, Computers and Mathematics with Applications 56 (2008) 1138–1145. [15] H. Wang, K. Wang, T. Sircar, A direct O(Nlog2N) finite difference method for fractional diffusion equations, Journal of Computational Physics 229 (2010) 8095–8104. [16] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, Journal of Computational Physics 205 (2005) 719–736. [17] M. Cui, Compact finite difference method for the fractional diffusion equation, Journal of Computational Physics 228 (2009) 7792–7804. [18] P. Zhuang, F. Liu, Implicit difference approximation for the time fractional diffusion equation, Journal of Applied Mathematics and Computing 22 (2006) 87–99. [19] F. Liu, S. Shen, V. Anh, I. Turner, Analysis of a discrete non-Markovian random walk approximation for the time fractional diffusion equation, ANZIAM Journal 46 (2005) 488–504. [20] C.-M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, Journal of Computational Physics 227 (2007) 886–897. [21] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, Journal of Computational Physics 225 (2007) 1533– 1552. [22] G.J. Fix, J.P. Roof, Least squares finite-element solution of a fractional order two-point boundary value problem, Computers Mathematics with Applications 48 (2004) 1017–1033. [23] Y. Jiang, J. Ma, High-order finite element methods for time-fractional partial differential equations, Journal of Computational and Applied Mathematics 235 (2011) 3285–3290. [24] C. Li, Z. Zhao, Y. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Computers and Mathematics with Applications 62 (2011) 855–875. [25] Q. Yang, I. Turner, F. Liu, M. Ilis, Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions, SIAM Journal on Scientific Computing 33 (2011) 1159–1180. [26] J.T. Katsikadelis, The BEM for numerical solution of partial fractional differential equations, Computers and Mathematics with Applications 62 (2011) 891–901. [27] H. Brunner, L. Ling, M. Yamamoto, Numerical simulations of 2D fractional subdiffusion problems, Journal of Computational Physics 229 (2010) 6613– 6622. [28] W. Chen, L. Ye, H. Sun, Fractional diffusion equations by the Kansa method, Computers and Mathematics with Applications 59 (2010) 1614–1620. [29] Q. Liu, Y. Gu, P. Zhuang, F. Liu, Y. Nie, An implicit RBF meshless approach for time fractional diffusion equations, Computational Mechanics 48 (2011) 1– 12. [30] Y. Gu, P. Zhuang, F. Liu, An advanced implicit meshless approach for the non-linear anomalous subdiffusion equation, CMES: Computer Modeling in Engineering and Sciences 56 (2010) 303–334. [31] S. Harald, Algorithm 368: numerical inversion of Laplace transforms [D5], Communications of the ACM 13 (1970) 47–49. [32] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Operations Research 14 (1966) 444–459. [33] W. Chen, Z.J. Fu, B.T. Jin, A truly boundary-only meshfree method for inhomogeneous problems based on recursive composite multiple reciprocity technique, Engineering Analysis with Boundary Elements 34 (2010) 196–205. [34] Z.J. Fu, W. Chen, W. Yang, Winkler plate bending problems by a truly boundary-only boundary particle method, Computational Mechanics 44 (2009) 757–763. [35] Z.J. Fu, W. Chen, C.Z. Zhang, Boundary particle method for Cauchy inhomogeneous potential problems, Inverse Problems in Science and Engineering 20 (2012) 189–207. [36] W. Chen, Z.J. Shen, L.J. Shen, G.W. Yuan, General solutions and fundamental solutions of varied orders to the vibrational thin, the Berger, and the Winkler plates, Engineering Analysis with Boundary Elements 29 (2005) 699–702. [37] W. Chen, Z.J. Fu, Q.H. Qin, Boundary particle method with high-order trefftz functions, CMC-Computers Materials and Continua 13 (2010) 201–217. [38] A.J. Nowak, A.C. Neves, The Multiple Reciprocity Boundary Element Method, Computational Mechanics Publication, 1994. [39] J. Abate, P.P. Valkó, Multi-precision Laplace transform inversion, International Journal for Numerical Methods in Engineering 60 (2004) 979–993. [40] J. Abate, W. Whitt, A unified framework for numerically inverting Laplace transforms, INFORMS Journal on Computing 18 (2006) 408–421. [41] X. Li, Convergence of the method of fundamental solutions for solving the boundary value problem of modified Helmholtz equation, Applied Mathematics and Computation 159 (2004) 113–125. [42] D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, New York, 1977. [43] M.A. Golberg, C.S. Chen, The method of fundamental solutions for potential, Helmholtz and diffusion problems. in: Boundary Integral Methods: Numerical and Mathematical Aspects, Computational Mechanics Publications, 1998, pp. 103–176. [44] A.H. Barnett, T. Betcke, Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains, Journal of Computational Physics 227 (2008) 7003–7026. [45] C.S. Chen, M.A. Golberg, Y.F. Rashed, A mesh free method for linear diffusion equations, Numerical Heat Transfer, Part B 33 (1998) 469–486.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Boundary particle method for Laplace transformed time fractional diffusion equations Zhuo-Jia Fu a,b, Wen Chen a,b,⇑, Hai-Tian Yang b a b
College of Mechanics and Materials and College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing, Jiangsu 210098, PR China State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China
a r t i c l e
i n f o
Article history: Received 19 March 2012 Received in revised form 1 October 2012 Accepted 18 October 2012 Available online 6 November 2012 Keywords: Boundary particle method Laplace transform Numerical inverse Laplace transform Meshless Time fractional derivative Anomalous diffusion
a b s t r a c t This paper develops a novel boundary meshless approach, Laplace transformed boundary particle method (LTBPM), for numerical modeling of time fractional diffusion equations. It implements Laplace transform technique to obtain the corresponding time-independent inhomogeneous equation in Laplace space and then employs a truly boundary-only meshless boundary particle method (BPM) to solve this Laplace-transformed problem. Unlike the other boundary discretization methods, the BPM does not require any inner nodes, since the recursive composite multiple reciprocity technique (RC-MRM) is used to convert the inhomogeneous problem into the higher-order homogeneous problem. Finally, the Stehfest numerical inverse Laplace transform (NILT) is implemented to retrieve the numerical solutions of time fractional diffusion equations from the corresponding BPM solutions. In comparison with finite difference discretization, the LTBPM introduces Laplace transform and Stehfest NILT algorithm to deal with time fractional derivative term, which evades costly convolution integral calculation in time fractional derivation approximation and avoids the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. Error analysis and numerical experiments demonstrate that the present LTBPM is highly accurate and computationally efficient for 2D and 3D time fractional diffusion equations. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Fractional diffusion equation [1] is considered a recent alternative model to describe anomalous diffusion phenomena [2] in a wide range of engineering and physics fields [3], such as electron transportation [4], seepage [5], magnetic plasma [6], dissipation [7] and turbulence [8], just to mention a few. Compared with the nonlinear models, the fractional anomalous diffusion models require fewer and easily available parameters to accurately describe the physical phenomena. Like the nonlinear models, the analytical solutions of fractional derivative models are largely not available, except for the cases with simple initial and boundary conditions by using Laplace transform [1,9] and Fourier–Laplace transforms. Hence the numerical simulation [10–12] plays an essential role in the analysis of fractional diffusion models. Nowadays the finite difference method (FDM) [13–15] is a popular and dominant numerical technique for time and spatial discretization of fractional derivative equations. Its convergence, accuracy, and stability have extensively been discussed in the literatures [16–19]. However, due to the memory (non-local) effect arising from the convolution integral expressions of time (space) fractional derivative term, the FDMs in these cases are no longer of a finite local discretization scheme and ⇑ Corresponding author at: College of Mechanics and Materials and College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing, Jiangsu 210098, PR China. Tel.: +86 25 83786873; fax: +86 25 83736860. E-mail address: [email protected] (W. Chen). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.10.018
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
53
encounter an exponential increase of computing costs with advancing time and large computational domain. Consequently, numerical results in most reports are restricted in 1D and 2D problems. This study considers time fractional derivative equations, which only has fractional derivative in time and integer differential operator in space. Most reports employ the finite difference method for temporal discretization but introduce the other numerical methods for spatial discretization, such as Fourier method [20], spectral method [21], finite element method [22–25], boundary element method [26], and radial basis function meshless collocation method [27–30]. In comparison with traditional FDMs, these methods can mitigate, to a certain extent, computing costs for large computational domain problems. But most of these methods still employ finite difference schemes for temporal discretization, and inevitably encounter an exponential increase of computing costs with advancing time and thus have low efficiency in the simulation of long timehistory fractional diffusion. In order to overcome the drawback of FDM for temporal discretization, this paper presents a novel approach, called Laplace transformed boundary particle method (LTBPM), for efficient modeling of time fractional diffusion equations. The LTBPM introduces Laplace transform to convert time fractional diffusion equation into a time-independent inhomogeneous equation in Laplace space. To retrieve the final time-dependent solutions from Laplace space solutions, the present LTBPM chooses a sophisticated Stehfest numerical inverse Laplace transform (NILT) [31], which is originated from Gaver [32] and applied to various scientific and engineering problems. It should be mentioned that, like the other NILT algorithms, the Stehfest algorithm is also an ill-posed problem and the magnification of truncation error results in a loss of numerical accuracy. For this reason, it is highly important to employ a highly-accurate and efficient spatial discretization method to solve the transformed time-independent diffusion equation in Laplace space. Fortunately, the authors recently develop boundary particle method (BPM) [33–35], which only requires boundary-only meshless discretization of inhomogeneous problem with integration-free merit and is very efficient and accurate for high-dimensional problems. In the present method, the BPM is used to solve the corresponding time-independent inhomogeneous equation in Laplace space. Therefore, the present LTBPM employs Laplace transform and Stehfest numerical inverse Laplace transform instead of the step-by-step FDMs for temporal discretization, which evades costly convolution integral calculation in time fractional derivative approximation and avoids the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. The rest of this paper is organized as follows. Section 2 introduces Laplace transformed boundary particle method (LTBPM) for numerical modeling of time fractional diffusion equations, followed by numerical results and discussions in Section 3. Some remarks are drawn upon the numerical results in Section 4.
2. Methodology 2.1. Time fractional diffusion problem Without loss of generality, we consider a time fractional diffusion equations in a bounded domain X with piecewise smooth boundary @ X ¼ CD þ CN ðCD \ CN ¼ ;Þ;
@ a uðx; tÞ ¼ DDuðx; tÞ þ ~ v ruðx; tÞ kuðx; tÞ þ Q ðx; tÞ; @t a
0 < a < 1;
x 2 X; t 2 ð0; TÞ;
ð1Þ
with boundary conditions
uðx; tÞ ¼ g 1 ðx; tÞ;
x 2 CD ;
@uðx; tÞ ¼ g 2 ðx; tÞ; @n
x 2 CN ;
t 2 ð0; TÞ;
ð2aÞ
t 2 ð0; TÞ;
ð2bÞ
and initial condition
uðx; 0Þ ¼ u0 ðxÞ;
x 2 X;
ð3Þ
where Q ðx; tÞ; g 1 ðx; tÞ; g 2 ðx; tÞ and u0 ðxÞ are known functions; x ¼ ðx; yÞ for 2D problem and x ¼ ðx; y; zÞ for 3D problem; D the diffusion coefficient, k the reaction coefficient, ~ v velocity vector, n the unit outward normal, T the total time to be considered, a and @t@ a the Caputo fractional derivative of order a with respect to t defined by
@ a uðx; tÞ 1 ¼ @t a Cð1 aÞ
Z 0
t
@uðx; gÞ dg ; @ g ðt gÞa
for more details see Podlubny [1].
0 < a < 1;
ð4Þ
54
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 1. The roadmaps of the LTBPM and the FDM for time fractional diffusion problems.
2.2. Implementation procedure Fig. 1 shows the roadmaps of the LTBPM and the FDM for time fractional diffusion problems. Unlike the one-step procedure in the FDM, the LTBPM includes three-step: (a) Laplace transform, (b) Boundary particle method (BPM), and (c) Numerical inverse Laplace transform (NILT). 2.2.1. Step (a) Laplace transform Laplace transform converts the time fractional diffusion problems from time domain to Laplace space domain. Let f(t) be a causal time domain function as an independent variable t P 0 and ~f ðsÞ denotes its image in Laplace space domain. Laplace transform is given by
~f ðsÞ ¼ Lðf ðtÞÞ ¼
Z
1
f ðtÞest dt;
ð5Þ
0
where s is Laplace transform parameter, and unless otherwise specified the quantities in Laplace space domain are denoted by an over tilde. Laplace transform of the Caputo fractional derivative [1] can be written as
a @ f ðtÞ ¼ sa ~f ðsÞ sa1 f ð0Þ; L a @t
0 < a < 1;
t P 0:
ð6Þ
Applying Laplace transform to Eqs. (1) and (2) produces
~ ðx; sÞ ¼ Fðx; sÞ; Ru ~ ðx; sÞ ¼ g~1 ðx; sÞ; u ~ ðx; sÞ @u ¼ g~2 ðx; sÞ; @n
0 < a < 1;
x 2 X;
x 2 CD ; x 2 CN ;
ð7Þ ð8aÞ ð8bÞ
where R ¼ DD þ ~ v r ðk þ sa Þ; Fðx; sÞ ¼ Q~ ðx; sÞ sa1 u0 ðxÞ. It is of interest to point out that, in comparison with finite difference discretization, Laplace transform evades costly convolution integral calculation of time fractional derivative and avoids the effect of time step on numerical accuracy and stability. 2.2.2. Step (b) Boundary particle method In the Laplace space, the boundary particle method (BPM) is implemented to solve Laplace-transformed time-indepen~ can be expressed as dent inhomogeneous problems (7) and (8) via meshless boundary nodes. The corresponding solution u
~p ; ~¼u ~h þ u u
ð9Þ
~ h and u ~ p are the homogeneous and the particular solutions, respectively. The particular solution u ~ p satisfies where u
~ p g ¼ F; Rfu
ð10Þ
~ h has to satisfy not only the but does not necessarily satisfy boundary conditions. In contrast, the homogeneous solution u corresponding homogeneous equation
~ h g ¼ 0; Rfu
ð11Þ
but also boundary conditions
~p ; ~ h jC ¼ g~1 u u D
ð12aÞ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~h ~ @u ¼ g~2 @ up : @n CN @n
55
ð12bÞ
~ h of Eqs. (11) and (12) can be efficiently calculated by boundary-type In the BPM procedure, The homogeneous solution u numerical techniques with the nonsingular general solutions [36], Trefftz functions [37] and singular fundamental solutions [38]. ~ p , the BPM employs the recursive composite multiple reciprocity method (RC-MRM) To evaluate the particular solution u [33] to annihilate the frequently encountered source term Fðx; sÞ in Eq. (10) by using composite differential operators
Lm . . . L2 L1 fFðx; sÞg ffi 0;
ð13Þ
where L1, L2, . . . Lm are differential operators of the same or different types, and Fðx; sÞ represents the polynomial, exponential and trigonometric functions or the combination of these functions in this study, which frequently appear in engineering problems. Consequently, Eq. (10) is equivalently transformed into the following high-order homogeneous problem
8 ~ p ðx; sÞ ¼ 0 x 2 X Lm . . . L2 L1 Ru > > > > .. < . ; > > L1 Ru ~ p ðx; sÞ ¼ L1 ðFðx; sÞÞ x 2 @ X > > : ~ p ðx; sÞ ¼ Fðx; sÞ Ru x 2 @X
0 < a < 1:
ð14Þ
Rearrange Eqs. (11) (12) and (14), we obtain the following high-order homogeneous problem
8 ~ ðx; sÞ ¼ 0 Lm . . . L2 L1 Ru > > > > > .. > > . > > > < ~ L1 Ruðx; sÞ ¼ L1 ðFðx; sÞÞ > ~ ðx; sÞ ¼ Fðx; sÞ > Ru > > > > ~ ðx; sÞ ¼ g~1 ðx; sÞ; > u > > > ~ ðx;sÞ : @u ¼ g~2 ðx; sÞ; @n
x2X
x 2 @X x 2 @X
0 < a < 1:
ð15Þ
x 2 CD x 2 CN
Then the solution of the transformed inhomogeneous Eqs. (7) and (8) is equal to that of the homogeneous problem (15). can be expressed as a linear combination of singular fundamental solutions with unHence, the approximate solution u known coefficients faij g as follows
ðxÞ ¼ u
m X N X aij Gi ðx yj Þ;
x 2 X;
ð16Þ
i¼0 j¼1
where N represents the number of source points, fyj g denotes the source points placed on the fictitious boundary outside the physical domain X, and Gi means the high-order singular fundamental solution of the differential operator Li (set L0 ¼ R). Some frequently-used high-order singular fundamental solutions are listed in Appendix. It is of interest to point out that the BPM is a truly boundary-only meshless and integration-free numerical technique. 2.2.3. Step (c) Numerical inverse Laplace transform Numerical inverse Laplace transform (NILT) is implemented to convert the BPM solutions from Laplace space domain to time domain. As the NILT is an ill-posed problem, the truncation error magnification appears and results in a loss of numerical accuracy. Here we employs the well-established Stehfest NILT algorithm originated from Gaver [32] to retrieve the timedependent solution of the original problem (1)–(3) from the BPM solutions in Laplace space. An approximate value fs of the Stehfest NILT of f(t) for a specific time t is given by
fs ðtÞ ¼
M ln 2 X ln 2 i ; V i ~f N t i¼1 t
ð17Þ
where ~f N is the Laplace space solutions, and M þi 2
V i ¼ ð1Þ
min ði;M 2Þ
X
k¼½iþ1 2
M
M 2
k 2 ð2kÞ! ; k !k!ðk 1Þ!ði kÞ!ð2k iÞ!
in which [C] denotes the nearest integers less than or equal to C. For a specific instant t, one should solve M boundary value problems (7) and (8) for the corresponding Laplace transform parameter s ¼ lnt 2 i; i ¼ 1; 2; . . . ; M. Then the solution of the time fractional diffusion Eqs. (1)–(3) can be retrieved by using the expression (17).
56
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
2.3. Error analysis This section presents error analysis of the present LTBPM. Since all of Laplace transforms are analytical in Step (a), and LTBPM’s numerical error arises largely in Step (b) and Step (c). First, we discuss the numerical error in Step (c), i.e. Stehfest numerical inverse Laplace transform (NILT). Abate et al. [39,40] did extensive numerical experiments to investigate the parameter effect on numerical accuracy, and summarized the following conclusion: ‘‘If p significant digits are desired, then let M be the positive integer [2.2p]; Given M, set the system precision at q ¼ ½1:1M; Given the transform ~f N and the argument t, then calculate fs by Eq. (17)’’. According to Abate’s conclusion, we obtain the following heuristic error estimation: ~
~
Remark 1. If the error of input data is 10ðqþ1Þ 6 kf N~f k 6 10q with even positive integer M via q ¼ ½1:1M, then the final kf k s f k error is 10ðpþ1Þ 6 kfkf 6 10p , where M ¼ ½2:2p. k It can be found from Remark 1 that the Stehfest algorithm tends to demand highly precise input data ~f N in order to yield good accuracy in the following NILT calculations. For example, if input data ~f N is single-precision data, namely, q = 7 with M = 8, then p = 4. If input data ~f N are double-precision data, namely, q = 16 with M = 16, then p = 7.In the present methodology, the input data ~f are obtained by solving time-independent inhomogeneous problems (7) and (8) in Laplace space. Therefore, it is necessary to employ a highly accurate numerical method in Step (b) to guarantee high accuracy in the final solutions. For this purpose, we implement the boundary particle method (BPM) with singular fundamental solutions in Step (b) to obtain highly accurate solution in Laplace space. As we will see in the following, the Laplace-transformed governing Eq. (7) degenerates to modified Helmholtz equation when the parameter ~ v ¼ ½0; 0. In BPM’s error analysis, we assume that the solution u~, the source term function F, and the ~1 are smooth enough on X. At first, the following homogeneous modified Helmholtz equation boundary condition function g with Dirichlet boundary condition under a unit disc is considered
~ ¼ 0; ðD k2 Þu
ð18Þ
~ jC ¼ g~1 : u D
ð19Þ
According to Theorem 2.3 in [41], we obtain a simplified error estimation N
~u kL1 ðXÞ 6 skg~1 kL1 ðC Þ ku D
d
N
1d
ð20Þ
;
is the corresponding approximation solution, the fictitious boundary parameter d > 1 represents the radius of a conwhere u centric disc Xd with unit disc X, N denotes the number of source points on @ Xd or collocation points on CD, s is a constant independent of N. The next step turns to consider the following inhomogeneous modified Helmholtz equation
~ ¼ F; ðD k2 Þu
ð21Þ
~ jC ¼ g~1 ; u D
ð22Þ
in which the source term F can be eliminated by using second-order differential operator L1, namely, L1 fFg ¼ 0 on X. ~ p and u ¼u h þ u p be the exact and numerical solution of the inhomogeneous problem (21) and (22), ~¼u ~h þ u Theorem 1. Let u Then p is the particular solution of Eq. (21), and ðD k2 Þu p ¼ F. where u N
~u kL1 ðXÞ 6 s1 kg~1 u p k L1 ðC Þ ku D where
d
N
1d
þ eF ;
ð23Þ
s1 is a constant independent of N, and eF satisfies
eF
8 > s2 kFkL1 ðCD Þ dN ; L1 ¼ D; > > < N ¼ s3 kFkL1 ðCD Þ d N ; L1 ¼ D k21 ; ; 1d > > > : s kFk 1 dN ; L1 ¼ D þ k21 4 L ðCD Þ
ð24Þ
in which s2, s3 and s4 are the constant independent of N. ^ h be the exact solution of the following problem Proof. Let u
^ h ¼ 0; ðD k2 Þu
ð25Þ
p : ^ h jC ¼ g~1 u u D
ð26Þ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
57
^¼u ^h þ u p ; then we have Let u
~u kL1 ðXÞ 6 ku ~u ^ kL1 ðXÞ þ ku ^u kL1 ðXÞ : ku
ð27Þ
By using Eq. (20), we have N
^u kL1 ðXÞ ¼ ku ^h u h kL1 ðXÞ 6 s1 kg~1 u p kL1 ðC Þ ku D
d
N
1d
ð28Þ
:
Notice that in X
~u ^Þ ¼ F F; ðD k2 Þðu
ð29Þ
and on CD
~u ^ÞjC ¼ g~1 ðg~1 u p þ u p Þ ¼ 0: ðu D
ð30Þ
According to a priori estimate [42], we obtain
1 : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ
ð31Þ
Consider the following problems
L1 v ¼ 0;
ð32Þ
v jC
ð33Þ
D
¼ F;
where F and F are its exact and numerical solution, respectively. For L1 ¼ D k21 , Eq. (31) can be rewritten by using Eq. (20)
1 6 s3 kFk 1 ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ
d
N N
1d
ð34aÞ
:
For L1 ¼ D, Eq. (31) can be rewritten by using Theorem 2 in [43]
1 6 s2 kFk 1 dN : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ For L1 ¼ D þ
k21 ,
ð34bÞ 2
Eq. (31) can be rewritten by using Theorem 3 with q > R in [44]
1 6 s4 kFk 1 dN : ~u ^ kL1 ðXÞ 6 skF Fk ku L ðXÞ L ðCD Þ
ð34cÞ
Therefore, the conclusion of the theorem follows from (27), (28) and (34). It can be observed from Eq. (23) that the BPM with singular fundamental solution has nearly exponential convergence for some inhomogeneous modified Helmholtz problems. Hence, the BPM error decreases exponentially as N ! 1 and d ! 1. As a price pays for such a high accuracy, the condition number of BPM interpolation matrix increases exponentially with N ! 1 and d ! 1. Fix on d, the BPM error tends to decrease exponentially and then oscillate slightly below a certain error level with an increase of N, due to the machine round-off error which becomes dominant in the final error at large value of N. Extensive numerical experiments [33,34,37] also demonstrate this observation. h 3. Numerical results and discussion In this section, the efficiency, accuracy and convergence of the proposed LTBPM are tested to four time fractional diffusion equations. The LTBPM solutions are compared with domain-type radial basis function (DRBF) solution [29] and analytical solution. For quantitative examinations, the notations of the maximum absolute and relative average error norms are defined below
ðjÞj; MerrðuÞ ¼ max juðjÞ u 16j6NT
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPNT u j¼1 ðuðjÞ u ðjÞÞ2 ; RerrðuÞ ¼ t PNT 2 j¼1 ðuðjÞÞ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPNT u j¼1 ðu;x ðjÞ u ;x ðjÞÞ2 ; Rerrðu;x Þ ¼ t PNT 2 j¼1 ðu;x ðjÞÞ
ð35aÞ
ð35bÞ
ð35cÞ
ðjÞ; u ;x ðjÞ and uðjÞ; u;x ðjÞ are the analytical and numerical solutions at xj , respectively, and NT denotes the total number where u of uniform test points in the domain of interest. In this study, NT is taken to be 121 in 2D cases and 1331 in 3D cases. In order to investigate the BPM accuracy in Step (b), we define
58
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~ Þ ¼ max ðMerrðu ~ ÞÞ; merrðu
ð36aÞ
~ ÞÞ; ~ Þ ¼ max ðRerrðu rerrðu
ð36bÞ
~;x ÞÞ; ~ ;x Þ ¼ max ðRerrðu rerrðu
ð36cÞ
16j6M
16j6M
16j6M
~; u ~ ;x are the solutions of Eqs. (7) and (8), i.e. the NILT input data ~f N in Eq. (17). where u Unless otherwise specified, the fictitious boundary parameter d = 6. For 2D cases, the source points are placed on a fictitious circle centering at geometric origin op and the fictitious boundary parameter d is its radius. For 3D cases, the fictitious boundary parameter d is defined as
d¼
yj xj : xj op
ð37Þ
Example 1. We first consider the time fractional advection-diffusion equation with the Dirichlet boundary conditions in a unit square domain X ¼ ½0; 1 ½0; 1. In order to compare with the reference results [29], the parameters in Eqs. (1)–(3) are h 2a i 2 xþy chosen as D ¼ 1; ~ v ¼ ½0; 0; k ¼ 0, Q ðx; tÞ ¼ C2tð3 ; g 1 ðx; tÞ ¼ t2 exþy ; u0 ðxÞ ¼ 0; where a ¼ 0:85; op = (0.5, 0.5) and aÞ 2t e the analytical solution uðx; tÞ ¼ t2 exþy , in which x ¼ ðx; yÞ 2 X. 2 exþy ; By applying the Laplace transform, we have the operators and functions R ¼ D sa ; Fðx; sÞ ¼ s43 þ s3 a g~1 ðx; sÞ ¼ s23 exþy in transformed problems (7) and (8). It is noted that the source term Fðx; sÞ can be annihilated by employing the differential operator L ¼ D 2 in the BPM.Table 1 displays the numerical accuracy of input data ~f obtained by the BPM N
1
with different boundary-only discretizations at T = 1. Numerical results verify that the BPM errors tend to decrease rapidly and then oscillate slightly below a certain error level with node refinements, which is in accordance with the BPM error analysis in Section 2.3. Fig. 2 shows the LTBPM error variations with respect to M in Eq. (17) by using different boundary node spacings (Dh ¼ 0:2; 0:1; 0:05) in comparison with the final errors obtained by the NILT with exact input data ~f . From Fig. 2, it can be found that all of the final errors tend to decrease rapidly and then increase with an increase of M, and all the turning points of error curves depend on the accuracy of NILT input data. It should be mentioned that in this case the NILT uses the BPM solutions ~f N (Dh ¼ 0:2; 0:1; 0:05) as input data to obtain the similar accuracy with exact input data ~f when M 6 14. It can be observed from Table 1 and Fig. 2 that when M 6 14 the loss of numerical accuracy mainly occurs in the NILT process rather than in the BPM computation. It is because the BPM solutions with Dh ¼ 0:2 are accurate enough to be the NILT input data, and the ill-posed NILT process results in the loss of numerical accuracy. This phenomenon can be also found in the method of particular solution with NILT for time-dependent problems in reference [45]. Furthermore, Table 2 presents the LTBPM errors with M ¼ 10 at T = 1 in comparison with the DRBF method. Numerical results show that the proposed LTBPM provides very accurate results by boundary-only discretization with node spacing Dh ¼ 0:2. In contrast, according to the convergence rate formula [29], the DRBF method requires whole domain discretization with smaller node spacing Dh ¼ 0:1 and very small time step Dt ¼ 0:004 to obtain the similar accuracy as the LTBPM. Namely, to achieve the similar accuracy, the LTBPM requires solving 20 20 linear systems with ðm þ 1Þ M ¼ 20 times, while the DRBF method requires solving 100 100 linear systems with T=Dt ¼ 250 times. Roughly speaking, the present study shows that the DRBF method requires more computational resources than the LTBPM. Moreover, in the time fractional derivative approximation, the FDM encounters an exponential increase of computing costs with advancing time because its memory requirements can grow in a dramatic way. In contrast, the LTBPM uses Laplace transform to evade this enormous effect on memory demand. These sharp comparisons indicate that the LTBPM is far more efficient for a long time-history fractional diffusion simulation. Example 2. The next case is a three-dimensional time fractional diffusion equation with the Dirichlet boundary conditions 2a
in a cube domain X ¼ ½0; 2 ½0; 2 ½0; 2. The parameters are Q ðx; tÞ ¼ ½C2tð3aÞ þ t 2 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ þ 6t2 , D ¼ 1; ~ v ¼ ½0; 0; k ¼ 1, u0 ðxÞ ¼ 0, g 1 ðx; tÞ ¼ t ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ, where a ¼ 0:7, op = (1, 1, 1) and the analytical solution is uðx; tÞ ¼ t2 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ, in which x ¼ ðx; y; zÞ 2 X. As in the first case, applying the Laplace transform yields the operators and source functions R ¼ D ðk þ sa Þ, 2 12 ~ 2 Fðx; sÞ ¼ ðs23 þ s3 1 ðx; sÞ ¼ s3 ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ in the transformed equaa Þ½xð2 xÞ þ yð2 yÞ þ zð2 zÞ þ s3 , g tions (7) and (8). The source term Fðx; sÞ can be annihilated by employing the differential operators L2 ¼ L1 ¼ D. 2
Table 1 The BPM errors with different boundary-only discretizations at T = 1. BPM
~Þ merrðu
~Þ rerrðu
~ ;x Þ rerrðu
Dh = 0.2 Dh = 0.1 Dh = 0.05
4.115e9 6.573e13 7.603e13
4.114e11 7.768e15 2.186e14
5.479e10 4.113e14 1.954e14
59
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 2. The final errors (a) Rerr(u) and (b) Rerr(u,x) with respect to M by using NILT with different input data (BPM solutions ~f N with Dh ¼ 0:2; 0:1; 0:05 and analytical solutions ~f ).
Table 2 The errors obtained by LTBPM with boundary-only discretization (Dh = 0.2) at T = 1 in comparison with the DRBF method with domain discretization (Dh = 0.1).
LTBPM DRBF(Dt = 0.1) [29] DRBF(Dt = 0.004)a a
Merr(u)
Rerr(u)
Rerr(u,x)
Rerr(u,y)
4.135e4 1.234e2 3.046e4
5.596e5 2.073e3 5.116e5
5.596e5 6.864e3 1.694e4
5.596e5 6.864e3 1.694e4
The DRBF results with Dt = 0.004 were obtained from the convergence rate formula OðDt2a Þ in reference [29].
Table 3 lists the BPM results at various time instants T by different boundary-only discretizations on the surface of 3D cube domain. It can be seen from Table 3 that the BPM with the same boundary discretization achieves the similar accuracy on various time instants T, even with small or large time instants T (T = 0.01, T = 100). It indicates that the present LTBPM avoids the effect of time step on numerical accuracy and stability encountered in the finite difference methods, and effectively simulates the long time-history fractional diffusion systems. Fig. 3 presents the LTBPM error variations with respect to M in Eq. (17) by using different boundary node spacings (Dh ¼ 0:4; 0:2) in comparison with the resulting errors by using exact input data in the NILT, and the corresponding LTBPM results at z = 1 cross section. This case verifies that the proposed LTBPM works equally well and efficiently for 3D fractional diffusion equations. It can be also observed from Fig. 3 that all of the final errors tend to decrease rapidly and then increase with increasing M, and all the turning points of error curves depend on the accuracy of NILT input data. Furthermore, the NILT uses the BPM solutions ~f N (Dh ¼ 0:2; 0:1; 0:05) as input data to obtain the similar accuracy with exact input data ~f when M 6 14. It also indicates that the final error is generated mainly by the Stehfest numerical inverse Laplace transform rather than by the BPM computation. It should be mentioned that Remark 1 is the heuristic error estimation. Hence we can estimate the optimal choice M is around 16 by using double-precision data type in the present numerical experiments. Moreover, the BPM numerical error in Step (b) (BPM NILT input data) is larger than the machine round-off error (exact NILT input data). Therefore, the optimal choice M may be less than 16 according to Remark 1. In addition, the numerical results in Examples 1 and 2 also demonstrate that the above discussions of the optimal choice M are basically correct, where M = 14 appears the best choice for simulation. So we choose M = 14 in the following numerical computation.
Table 3 The BPM errors with different boundary-only discretizations at various time instants T. BPM
Dh = 0.4
Dh = 0.2
T = 0.01 T=1 T = 100 T = 0.01 T=1 T = 100
~Þ rerrðu
~ ;x Þ rerrðu
1.045e9 1.045e9 1.044e9 4.475e12 4.455e12 2.967e12
7.464e9 7.464e9 7.457e9 1.986e11 2.666e11 3.147e11
60
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 3. LTBPM error variations with respect to M on different time instants (a) T = 0.01, (c) T = 1, (e) T = 100; and LTBPM solution at z = 1 cross section on different time instants (b) T = 0.01, (d) T = 1, (f) T = 100.
Example 3. This case tests the proposed LTBPM approach to the time fractional diffusion equation with zero Dirichlet/mixed boundary conditions in a square domain X ¼ ½1; 1 ½1; 1 [27]. The parameters in Eqs. (1)–(3) are chosen as D ¼ 1; ~ v ¼ ½0; 0; k ¼ 0, Q ðx; tÞ ¼ 0; g 1 ðx; tÞ ¼ g 2 ðx; tÞ ¼ 0, op = (0, 0), in which x ¼ ðx; yÞ 2 X. For Dirichlet boundary conditions, the analytical solution associated with initial condition
u0 ðxÞ ¼ cos is given by
p p x cos y ; 2 2
ð38Þ
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
61
Fig. 4. The maximum norm of the analytical solution with Dirichlet and mixed boundary conditions, respectively, for a ¼ 0:1; 0:2; . . . ; 0:9 within time period T 2 ½0; 1.
Fig. 5. The final errors obtained by the NILT with exact input data for Dirichlet and mixed boundary conditions, respectively, at a ¼ 0:1; 0:2; . . . ; 0:9 within time period T 2 ½0; 1.
p p 1 uðx; tÞ ¼ Ea p2 t a cos x cos y : 2 2 2
ð39Þ
For mixed boundary conditions having Dirichlet boundary on x ¼ 1 and x ¼ 1 and Neumann boundary on y ¼ 1 and y ¼ 1. The initial condition is
u0 ðxÞ ¼ cos
p p x sin y : 2 2
ð40Þ
Then we have the analytical solution
p p 1 uðx; tÞ ¼ Ea p2 t a cos x sin y : 2 2 2
ð41Þ
where the one-parameter Mittag–Leffler function is defined as
Ea ðzÞ ¼
1 X
zj
j¼0
Cðaj þ 1Þ
; a > 0:
ð42Þ
62
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
~ Þ) with Fig. 6. LTBPM error profiles (Rerr(u)) with boundary-only discretizations (a) Dh = 0.2, (c) Dh = 0.1; and the corresponding BPM error profiles (rerrðu boundary-only discretizations (b) Dh = 0.2, (d) Dh = 0.1 for Dirichlet problem in Example 3.
As in the previous cases, we employ the Laplace transform to get the expression of R ¼ D sa ; Fðx; sÞ ¼ sa1 u0 ðxÞ; ~2 ðx; sÞ ¼ 0 in transformed Eqs. (7) and (8). The source term Fðx; sÞ can be annihilated by employing the differential g~1 ðx; sÞ ¼ g 2 operator L1 ¼ D þ p2 in the BPM. It should be stressed that the BPM must add a constraint condition to guarantee the uniqueness of the solution of the transformed high-order homogeneous equations. The simplest way is to add the approximate formulation (16) to satisfy the governing Eq. (7) in few additional inner points. In this case, only one additional inner point (0, 0) is enough to get accurate and stable solutions. Fig. 4 displays the maximum norm of the analytical solutions with Dirichlet and mixed boundary condition problems, respectively, for a ¼ 0:1; 0:2; ; 0:9 and time period T 2 ½0; 1. It can be found that the maximum norms drop very quickly at the beginning instants, especially for small a. As time advances, the solution varies far less after a short period. Fig. 5 shows the final errors obtained by the NILT with exact input data for Dirichlet and mixed boundary condition problems, respectively, at a ¼ 0:1; 0:2; ; 0:9 and time period T 2 ½0; 1. We can see from Fig. 5 that the errors decrease with an decrease of fractional order a. ~ Þ) by using boundary-only Fig. 6 shows the LTBPM error profiles (RerrðuÞ) and the corresponding BPM error profiles (rerrðu discretization with node spacing Dh ¼ 0:2 and Dh ¼ 0:1 for Dirichlet boundary condition case under a ¼ 0:1; 0:2; ; 0:9 within T 2 ½0; 1. From Fig. 6, one can observe that the LTBPM achieves acceptable accuracy merely with Dh ¼ 0:2 boundaryonly discretization. With a further node refinement (Dh ¼ 0:1), the LTBPM almost obtains the same accuracy as the NILT with exact input data ~f as shown in Fig. 5. This means that the BPM solutions ~f N with Dh ¼ 0:1 are accurate enough to be the NILT input data, and the loss of numerical accuracy is mainly generated in the ill-posed NILT process due to its error magnification.
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
63
~ Þ) with Fig. 7. LTBPM error profiles (Rerr(u)) with boundary-only discretizations (a) Dh = 0.2, (c) Dh = 0.1; and the corresponding BPM error profiles (rerrðu boundary-only discretizations (b) Dh = 0.2, (d) Dh = 0.1 for mixed boundary problem in Example 3.
~ Þ) by using Fig. 7 further illustrates the LTBPM error profiles (RerrðuÞ) and the corresponding BPM error profiles (rerrðu boundary-only discretization with node spacing Dh ¼ 0:2 and Dh ¼ 0:1 for mixed boundary condition case under a ¼ 0:1; 0:2; . . . ; 0:9 within T 2 ½0; 1. The similar observations can be made as in the previous Dirichlet problem. Example 4. The last case is concerned with the fractional subdiffusion–convection problem with mixed boundary conditions v ¼ ½0:005; 0; k ¼ 0, Q ðx; tÞ ¼ 0, in a square domain X ¼ ½1; 1 ½1; 1. The corresponding parameters are D ¼ 1; ~ 5
u0 ðxÞ ¼ ðxþ1Þ , op = (0, 0), The Dirichlet boundary conditions g 1 ðx; tÞ ¼ x þ 1 on x ¼ 1 and x ¼ 1 and the Neumann boundary 16 conditions g 2 ðx; tÞ ¼ 0 on y ¼ 1 and y ¼ 1, in which x ¼ ðx; yÞ 2 X. @ We can derive the operators and source functions of the Laplace transformed Eqs. (7) and (8), R ¼ D þ 0:005 @x sa ; ~ Fðx; sÞ ¼ sa1 u0 ðxÞ; g~1 ðx; sÞ ¼ xþ1 ; g ðx; sÞ ¼ 0. It is easily to find that the source term Fðx; sÞ can be annihilated by 2 s employing the differential operators L3 ¼ L2 ¼ L1 ¼ D in the BPM. Due to the symmetry of this problem, Fig. 8 shows LTBPM solutions with boundary-only discretization (Dh ¼ 0:2) at the cross section (parallel to the x1-axis) with varied fractional order a for every 1/32 stepping up from the lower-right towards the diagonal. The proposed LTBPM utilizes d ¼ 2:5 and coarse discretization (Dh ¼ 0:2) to produce the similar solution profiles as given in the reference method [27]. It can be observed from Fig. 8 that the effect of convection is significant in the case of a ¼ 0:9, and the inflection points can more clearly be seen at the initial instants. As the fractional order a decreases,
64
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Fig. 8. LTBPM solution (cross section) of the fractional subdiffusion-convection problem with varied order a of fractional derivative by using boundary-only discretization (Dh = 0.2).
the effect of convection is weakened. It should be noting that the numerical solutions for different a vary tremendously near the left endpoint x ¼ 1. The above numerical results favourably suggest that sensors should be placed somewhere in 1 < x < 0 instead of in 0 < x < 1 in the experiment design to capture the dramatic variance of the solution.
4. Conclusions This paper develops a Laplace transformed boundary particle method for numerical modeling of time fractional diffusion equations. In contrast to the standard FDMs, the proposed scheme implements Laplace transform and Stehfest numerical Laplace inversion to handle the time fractional derivative. It first uses Laplace transform to convert time fractional diffusion equation into a time-independent inhomogeneous equation in Laplace space. And then it uses a truly boundary-only meshless BPM to discretize the spatial derivatives of this Laplace-transformed inhomogeneous equation only with boundary discretization nodes. Finally, it chooses Stehfest numerical Laplace inversion to retrieve the numerical solutions of time fractional diffusion equations from the corresponding BPM solutions in Laplace space. In comparison with the finite difference discretization, the present LTBPM introduces Laplace transform technique to evade costly convolution integral calculation in time fractional derivation approximation and to avoid the effect of time step on numerical accuracy and stability. Consequently, it can effectively simulate long time-history fractional diffusion systems. On the other hand, the BPM makes the present strategy very attractive to high-dimensional problems. Error analysis and numerical experiments demonstrate that the LTBPM is a competitive numerical method in the solution of 2D and 3D time fractional diffusion equations with the merits of mathematically simple, easy-to-program, inherently meshless, highly accurate, boundary-only discretization and integration-free.
65
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
Moreover, it is observed from numerical results that the Stehfest numerical inverse Laplace transform has the error magnification in ill-posed process results, which results in loss of numerical accuracy. To solve this ill-posed issue, the present LTBPM implements the BPM to produce the highly accurate NILT input data and guarantee that the NILT still achieves high accuracy in the final results. Furthermore, it is of interest to point out that the multiple precision arithmetic [39] becomes a popular and active research area in science and engineering community. It can be introduced to further overcome the loss of numerical accuracy in the ill-posed NILT process, which is now under intense study. The further results will be reported in a subsequent paper. Acknowledgements The authors thank Prof. CS Chen and the anonymous reviewers of this article for their very helpful comments and suggestions to significantly improve the academic quality of this article. This work was supported by National Basic Research Program of China (973 Project No. 2010CB832702), R&D Special Fund for Public Welfare Industry (Hydrodynamics, Grant No. 201101014), National Science Funds for Distinguished Young Scholars (Grant No. 11125208) and Programme of Introducing Talents of Discipline to Universities (111 project, Grant No. B12032). The first author would like to thank Hohai University Training Program for Excellent Doctoral Dissertation (Grant No. 2011B14814), Jiangsu Province Graduate Students Research and Innovation Plan (Grant No. CX10B_203Z) for financial support. Appendix A In the appendix, D denotes Laplacian operator, r the gradient operator, D the diffusivity coefficient, k a real number qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 j~ vj 2 l ¼ ð2D Þ þ kD , j foundation stiffness, and r
v and ~r the velocity vector and distance vector, known as the wave number, ~
the Euclidean norm between the point x and the origin, i.e. r ¼ kxk. Ym and Km represent the m-order Bessel and modified Bessel functions of the second kind, respectively. Note that the fundamental solution to a differential operator is not unique. For the Laplacian, a constant may be included. However, it is often omitted for it does not have a significant effect on the accuracy of the numerical results. As discussed in [38], the high-order fundamental solution of the Laplacian operator
Dmþ1 can be derived as follow ( r2m ðC m ln r Bm Þ; 2p /m ðxÞ ¼ F 1 r 2m1 ; ð2mÞ! 4p
x 2 R2 x 2 R3
ðA1Þ
;
Cm 1 where C 0 ¼ 1; B0 ¼ 0; C mþ1 ¼ 4ðmþ1Þ 2 ; Bmþ1 ¼ 4ðmþ1Þ2
Cm mþ1
þ Bm .
The high-order fundamental solution of Helmholtz-type operator ðD þ k2 Þmþ1 [33] is given by mþ1n=2 /m Y m1þn=2 ðkrÞ; F ðxÞ ¼ Am ðkrÞ
x 2 Rn :
ðA2Þ 2 mþ1
The high-order fundamental solution of modified Helmholtz-type operator ðD k Þ
/m F ðxÞ ¼ Am ðkrÞ
mþ1n=2
K m1þn=2 ðkrÞ;
[33] is given by
x 2 Rn :
ðA3Þ 2 mþ1
The high-order fundamental solution of modified convection-diffusion-type operator ðDD þ v r k Þ
[33] is given
by ~ v~r
mþ1n=2 /m e 2D K m1þn=2 ðlrÞ; x 2 Rn ; F ðxÞ ¼ Am ðlrÞ
ðA4Þ
2
where Am ¼ Am1 =2mk ; A0 ¼ 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Reports 339 (2000) 1–77. R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000. H. Scher, E.W. Montroll, Anomalous transit-time dispersion in amorphous solids, Physical Review B 12 (1975) 2455–2477. R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini, P. Paradisi, Discrete random walk models for space-time fractional diffusion, Chemical Physics 284 (2007) 521–541. D. del-Castillo-Negrete, B.A. Carreras, V.E. Lynch, Front dynamics in reaction-diffusion systems with levy flights: a fractional diffusion approach, Physical Review Letters 91 (2003) 018302. T.L. Szabo, J. Wu, A model for longitudinal and shear wave propagation in viscoelastic media, Journal of Acoustical Society of America 107 (2000) 2437– 2446. I.M. Sokolov, J. Klafter, A. Blumen, Ballistic versus diffusive pair dispersion in the Richardson regime, Physical Review E 61 (2000) 2717–2722. R. Gorenflo, Y. Luchko, F. Mainardi, Wright functions as scale-invariant solutions of the diffusion-wave equation, Journal of Computational and Applied Mathematics 118 (2000) 175–191. K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Computer Methods in Applied Mechanics and Engineering 194 (2005) 743–773. S. Momani, An algorithm for solving the fractional convection–diffusion equation with nonlinear source term, Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1283–1290.
66
Z.-J. Fu et al. / Journal of Computational Physics 235 (2013) 52–66
[12] P.P. Valko, J. Abate, Numerical inversion of 2-D Laplace transforms applied to fractional diffusion equations, Applied Numerical Mathematics 53 (2005) 73–88. [13] C. Tadjeran, M.M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, Journal of Computational Physics 220 (2007) 813–823. [14] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equations, Computers and Mathematics with Applications 56 (2008) 1138–1145. [15] H. Wang, K. Wang, T. Sircar, A direct O(Nlog2N) finite difference method for fractional diffusion equations, Journal of Computational Physics 229 (2010) 8095–8104. [16] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, Journal of Computational Physics 205 (2005) 719–736. [17] M. Cui, Compact finite difference method for the fractional diffusion equation, Journal of Computational Physics 228 (2009) 7792–7804. [18] P. Zhuang, F. Liu, Implicit difference approximation for the time fractional diffusion equation, Journal of Applied Mathematics and Computing 22 (2006) 87–99. [19] F. Liu, S. Shen, V. Anh, I. Turner, Analysis of a discrete non-Markovian random walk approximation for the time fractional diffusion equation, ANZIAM Journal 46 (2005) 488–504. [20] C.-M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, Journal of Computational Physics 227 (2007) 886–897. [21] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, Journal of Computational Physics 225 (2007) 1533– 1552. [22] G.J. Fix, J.P. Roof, Least squares finite-element solution of a fractional order two-point boundary value problem, Computers Mathematics with Applications 48 (2004) 1017–1033. [23] Y. Jiang, J. Ma, High-order finite element methods for time-fractional partial differential equations, Journal of Computational and Applied Mathematics 235 (2011) 3285–3290. [24] C. Li, Z. Zhao, Y. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Computers and Mathematics with Applications 62 (2011) 855–875. [25] Q. Yang, I. Turner, F. Liu, M. Ilis, Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions, SIAM Journal on Scientific Computing 33 (2011) 1159–1180. [26] J.T. Katsikadelis, The BEM for numerical solution of partial fractional differential equations, Computers and Mathematics with Applications 62 (2011) 891–901. [27] H. Brunner, L. Ling, M. Yamamoto, Numerical simulations of 2D fractional subdiffusion problems, Journal of Computational Physics 229 (2010) 6613– 6622. [28] W. Chen, L. Ye, H. Sun, Fractional diffusion equations by the Kansa method, Computers and Mathematics with Applications 59 (2010) 1614–1620. [29] Q. Liu, Y. Gu, P. Zhuang, F. Liu, Y. Nie, An implicit RBF meshless approach for time fractional diffusion equations, Computational Mechanics 48 (2011) 1– 12. [30] Y. Gu, P. Zhuang, F. Liu, An advanced implicit meshless approach for the non-linear anomalous subdiffusion equation, CMES: Computer Modeling in Engineering and Sciences 56 (2010) 303–334. [31] S. Harald, Algorithm 368: numerical inversion of Laplace transforms [D5], Communications of the ACM 13 (1970) 47–49. [32] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Operations Research 14 (1966) 444–459. [33] W. Chen, Z.J. Fu, B.T. Jin, A truly boundary-only meshfree method for inhomogeneous problems based on recursive composite multiple reciprocity technique, Engineering Analysis with Boundary Elements 34 (2010) 196–205. [34] Z.J. Fu, W. Chen, W. Yang, Winkler plate bending problems by a truly boundary-only boundary particle method, Computational Mechanics 44 (2009) 757–763. [35] Z.J. Fu, W. Chen, C.Z. Zhang, Boundary particle method for Cauchy inhomogeneous potential problems, Inverse Problems in Science and Engineering 20 (2012) 189–207. [36] W. Chen, Z.J. Shen, L.J. Shen, G.W. Yuan, General solutions and fundamental solutions of varied orders to the vibrational thin, the Berger, and the Winkler plates, Engineering Analysis with Boundary Elements 29 (2005) 699–702. [37] W. Chen, Z.J. Fu, Q.H. Qin, Boundary particle method with high-order trefftz functions, CMC-Computers Materials and Continua 13 (2010) 201–217. [38] A.J. Nowak, A.C. Neves, The Multiple Reciprocity Boundary Element Method, Computational Mechanics Publication, 1994. [39] J. Abate, P.P. Valkó, Multi-precision Laplace transform inversion, International Journal for Numerical Methods in Engineering 60 (2004) 979–993. [40] J. Abate, W. Whitt, A unified framework for numerically inverting Laplace transforms, INFORMS Journal on Computing 18 (2006) 408–421. [41] X. Li, Convergence of the method of fundamental solutions for solving the boundary value problem of modified Helmholtz equation, Applied Mathematics and Computation 159 (2004) 113–125. [42] D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, New York, 1977. [43] M.A. Golberg, C.S. Chen, The method of fundamental solutions for potential, Helmholtz and diffusion problems. in: Boundary Integral Methods: Numerical and Mathematical Aspects, Computational Mechanics Publications, 1998, pp. 103–176. [44] A.H. Barnett, T. Betcke, Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains, Journal of Computational Physics 227 (2008) 7003–7026. [45] C.S. Chen, M.A. Golberg, Y.F. Rashed, A mesh free method for linear diffusion equations, Numerical Heat Transfer, Part B 33 (1998) 469–486.
