An implicit MLS meshless method for 2-D time dependent fractional diffusion–wave equation

Accepted Manuscript An implicit MLS meshless method for 2-D time dependent fractional diffusionwave equation J.Y. Yang, Y.M. Zhao, N. Liu, W.P. Bu, T.L. Xu, Y.F. Tang PII: DOI: Reference:

S0307-904X(14)00395-3 http://dx.doi.org/10.1016/j.apm.2014.08.005 APM 10109

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

16 October 2013 24 March 2014 14 August 2014

Please cite this article as: J.Y. Yang, Y.M. Zhao, N. Liu, W.P. Bu, T.L. Xu, Y.F. Tang, An implicit MLS meshless method for 2-D time dependent fractional diffusion-wave equation, Appl. Math. Modelling (2014), doi: http:// dx.doi.org/10.1016/j.apm.2014.08.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An implicit MLS meshless method for 2-D time dependent fractional diffusion-wave equation JY Yanga , YM Zhaob, N Liua , WP Bua , TL Xuc , YF Tanga,∗ a LSEC,

ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China b School of Mathematics and Statistics, Xuchang University, Henan 461000, China c Department of Mathematics, School of Information,Renmin University of China, Beijing 100872, China

Abstract It is well accepted that fractional partial differential equations (FPDE) can be used to model many processes for which the normal partial differential equations (PDE) fail to describe precisely. Numerical approaches seem to be promising alternatives when exact solution of FPDE is difficult to derive. However, numerical solution of FPDE encounters new challenges brought in by the fractional order derivatives. In this paper we consider the 2D time dependent fractional diffusion-wave equation (FDWE) with Caputo derivative in temporal direction. We discretize the fractional order derivative with finite difference method (FDM) and present a moving least squares (MLS) meshless approximation in spatial directions which can be used to handle more complex problem domain. The convergence and stability properties of semi-discretized scheme related to time are theoretically analyzed. Finally, we conduct several numerical experiments to test our method for both regular and irregular node point distribution on rectangular and circular domain. The results indicate that the proposed method is accurate and efficient. Keywords: Fractional diffusion-wave equation, Caputo derivative, Finite difference method, MLS method

1. Introduction In the past decades fractional calculus has attracted considerable interests in the fields of engineering and science. One can refer to the treatises [1–3] for elementary knowledge of fractional calculus. Roughly speaking fractional space derivatives usually can be used to model anomalous diffusion or dispersion while fractional time derivatives to model some processes with ‘memory’ effects. As a counterpart of traditional integer order differential equation, fractional differential equation can be obtained by replacing the integer order derivatives with fractional ones in integer order differential equation. Fractional differential equations such as fractional kinetic equations, fractional diffusion equation, fractional advection-diffusion equation and fractional Foker-Planck equation can be used to model the transport dynamics in complex systems. Aleroev et al. [4] presented a new initial value problem of a nonlinear fractional differential equation with variable coefficient to model the seepage of a liquid to a chink in the cracked deformable layer. Specifically when replacing the first or second order time derivative in diffusion or wave equation respectively with fractional (α > 0 order) time derivative we get time fractional diffusion equation as 0 < α < 1 and time FDWE as 1 < α < 2. Specially, time FDWE can be viewed as the interpolation of diffusion and wave equation and be used to model the process which exhibits heat diffusion and wave propagation properties. Nigmatullin [5, 6] reported that time FDWE can be used to accurately model many of the universal electromagnetic, acoustic, and mechanical responses. Many theoretical work has been conducted on the time fractional diffusion equation and time FDWE. Wess and Schneider [7, 8] presented solution of FDWE in terms of Fox’s H-function. Mainardi [9, 10] obtained the fundamental solution of FDWE using Laplace transform in terms of auxiliary function M(z, β) with z = |x|/tβ denoting the similarity variables. Agrawal [11, 12] derived a general solution to FDWE with fourth order space derivative in bounded domain. ∗ Corresponding

author. Email address: [email protected] (YF Tang)

Preprint submitted to Applied Mathematical Modelling

August 23, 2014

In [13] Agrawal reduced the FDWE into a set of equations using the method of separation of variables and then used Laplace transform technique to obtain the fractional Green’s function and to describe the system’s response of a Duhamel integral type. Metzler and Klafler [14] solved the fractional diffusion equation with absorbing and reflecting boundary conditions using Fourier-Laplace transform and the the method of separation of variables. So far many researchers devoted themselves to numerical treatment of fractional diffusion equation and time FDWE. Yuste and Acedo [15, 16] presented a forward Euler scheme and weighted average finite difference scheme for time fractional diffusion equation and analyzed the stability of their method. However, they did not derive the convergence order of the proposed method. Liu et al. [17] considered the same equation with a first order discretization of time and space derivative and proved the stability. Zhuang and Liu [18] presented a O(τ + h2 ) scheme (τ and h denote stepsize in time and space direction respectively here and in the remainder of this section) for the fractional sub-diffusion equation and gave convergence and stability analysis. Xu et al. [19, 20] proposed spectral method for fractional sub-diffusion equation. In [19] they approximated the time derivative with finite difference method while space derivative with Legendre spectral method and proved the convergence and stability of their scheme. In [20] the authors discretized both time and space derivatives with spectral method and got a scheme of 2 − α (α represents the order of time derivative in the equation considered) order in time direction and with spectral accuracy in space direction. However, although with spectral accuracy, the spectral method suffers from a drawback that it imposes so many requirements on the smoothness of the exact solution. Sun [21] presented a O(τ3−α + h2 ) full discrete difference approximation for 1-D time FDWE by introducing two new intermediate variables and proved the convergence and stability of the proposed method. Tadjeran et al. [22] considered the fractional diffusion equation with RiemmanLiouville derivative in time direction and constructed a second order scheme based on Crank-Nicholson method and Richardson extrapolation. Unfortunately neither FDM nor spectral method could handel complex problem domain easily. The meshless methods are developed to overcome this drawback and become increasingly popular in the field of computational mechanics. Unlike FDM, finite element method (FEM) and spectral method, the meshless methods use a set of quasirandom node points to discretize the problem domain rather than grid or mesh structure demanded by FDM, FEM and spectral methods. This helps to cut down the computation cost in grid generation which is a crucial step in the implementation of FDM and FEM. Furthermore the quasi-random distribution of node points makes the meshless methods prone to design of some adaptive algorithms. Another distinguished advantage of the meshless methods over FDM and FEM is the accuracy due to the use of higher order meshless shape function. Generally the meshless methods fall into two categories. The first category is collocation based methods which include smooth particle hydrodymacis (SPH) [23], the moving least squares method [24], meshless collocation method [25] and finite point method (FPM) [26] etc. This kind of methods are introduced and developed relatively early and use Dirac-delta test function. The second category uses various global weak-forms and local weak-forms and its birth can date back to 20 years ago. This kind of meshless methods include the reproducing kernel particle method, the element-free Galerkin (EFG) method, the point interpolation method (PIM), the diffuse element method (DEM), the meshless local PetrovGalerking (MLPG) method, the local radial point interpolation method (LRPIM), the boundary node method (BNM) and the boundary point interpolation method (BPIM) (see [24] and references therein for the detailed discription of these methods). Gu [27] made a review on several kinds of meshless methods and compared them in terms of accuracy, convergence and efficiency. Application of the meshless methods in numerical simulation of FPDE is still in its first stage of development. Gu et al. [28] presented an implicit meshless approach based on the radial basis function (RBF) for the 2-D anomalous subdiffusion equation and discussed the convergence and stability of the space semi-discretized scheme. Zhuang et al. [29] constructed an implicit MLS meshless approximation for the 2-D fractional advection-diffusion equation with the order of time Caputo derivative in the interval of (0, 1) and analyzed the convergence and stability of the proposed scheme as well. Liu et al. [30, 31] solved the time fractional diffusion equation and fractal mobile/immobile transport model with RBF meshless approach and proved the convergence and stability of their proposed methods theoretically. To the extent of our knowledge, numerical treatment of the time FDWE with time Caputo derivative using the meshless methods has not been reported yet. It is well known that when 1 < α < 2 the Caputo derivative contains a second order derivative in the integrand which makes FDWE more difficult to handle than the fractional diffusion equation. This is the motivation of our present work. In this article we propose an implicit MLS meshless method for the 2-D time FDWE on both rectangular and circular domain. We discretize the time Caputo derivative with finite difference method combined with the MLS meshless method in space direction. The convergence and stability of our newly 2

developed approach is proved theoretically. To validate our method we conduct several numerical experiments. The remainder of this paper is organized as follows. In section 2 we present the semi-discretization in time direction and prove the convergence and stability of our scheme. In section 3 the MLS meshless approach is introduced and then used to derive a full discretization of the time-dependent FDWE. In section 4 we validate the scheme by several numerical examples on both rectangular and circular domain. In the final section we draw some conclusions. 2. Semi-discretization and analysis of convergence and stability Consider the two dimensional fractional diffusion-wave equation ∂α u(x, t) = c4u(x, t) + f (x, t), x ∈ Ω ⊂ R2 , 0 < t ≤ T, ∂tα

(1)

together with the initial conditions u(x, 0) = ϕ(x),

∂u(x, 0) = ψ(x), x ∈ Ω, ∂t

and the boundary conditions u(x, t) = g(x, t), x ∈ ∂Ω, t > 0, where c is the diffusion coefficient and f (x, t) denotes the source term, ∆ is the Laplacian operator, x = (x, y)T and t are space and time variables and Z t 2 1 ∂α u 1−α ∂ u(x, s) = (t − s) ds, (1 < α < 2), (2) ∂tα Γ(2 − α) 0 ∂s2 is Caputo derivative with respect to t with Γ being the gamma function. Rewriting the Caputo derivative gives [32] Z t 1 ∂α u ∂ ∂u(x, s) ∂α−1 ∂ = (t − s)−(α−1) ds = α−1 u(x, t). (3) α ∂t Γ(1 − (α − 1)) 0 ∂s ∂s ∂t ∂t Based on formula (3), the following theorem holds true. Theorem 1. Equation (1) is equivalent to the following partial integro-differential equation " 2 # Z t ∂ u(x, s) ∂2 u(x, s) c (t − s)α−2 + ds + F(x, t), ut (x, t) = ψ(x) + Γ(α − 1) 0 ∂x2 ∂y2

(4)

β

where F(x, t) = Jtα−1 f (x, t), and Jt is the βth order Riemann-Liouville fractional integral operator defined as Z t 1 (t − s)β−1 f (x, s)ds. Jtβ f (x, t) = Γ(β) 0

(5)

Proof. The statement can be easily proved by acting integral operator Jtα−1 on both sides of equation (1). Here we will not cover the detail. Now we settle down to derive a finite difference semi-discretization for equation (4). Let tn = nτ, where τ = T/N is the time stepsize. We also need the following lemma.

3

Lemma 1. [33]. For f (t) = tβ−1 g(t), with β > 0, and g(t) a smooth function, we have Jtα−1 f (tn ) = τα−1

n X

ωk f (tn−k ) + O(τ),

(6)

k=0

where ωk are the convolution weights determined by the generation function σ(ζ) = (1 − ζ)1−α . Replacing the time derivative in (4) with first order difference and using Lemma 1 we get u(x, tn+1 ) − u(x, tn ) = cτα

n+1 X

ωk ∆u(x, tn+1−k ) + τψ(x) + τF(x, tn+1 ) + Rn+1 ,

(7)

k=0

where Rn+1 = O(τ2 ). Let un denote the approximation of u(x, tn ). Then the following semi-discretization of (4) can be readily obtained un+1 − un = cτα

n+1 X

ωk ∆un+1−k + τψ + τF n+1 .

(8)

k=0

In order to study the convergence and stability of scheme (8) we introduce the inner product " (u, v) = u(x)v(x)dxdy, Ω

and L2 norm kvk2 = (v, v)1/2 =

""

v(x)2 dxdy Ω

#1/2

.

We also need the following two lemmas. Lemma 2. [32] The convolution weights ωn , n = 0, 1, 2, · · · , in (8) satisfy the following properties ωn ≥ 0, ωn − ωn+1 ≥ 0, ωn+2 − 2ωn+1 + ωn ≥ 0, n = 0, 1, 2, · · · . Lemma 3. [34] Let {bn }∞ n=0 be a sequence of real numbers which satisfy the following properties bn ≥ 0,

bn+1 − bn ≤ 0,

bn+1 − 2bn + bn−1 ≥ 0.

Then for any positive integer J and real vector (V1 , V2 , · · · , V J ) with J real entries,   n J−1 X X    b p Vn+1−p  Vn+1 ≥ 0. n=0

p=0

Based on Lemma 2 and 3 we can readily get the results of convergence and stability of scheme (8). Theorem 2. Assume that utt is bounded on the domain Ω × [0, T ], and u(x, t) and un are the solutions of equation (4) and semi-discretization scheme (8) respectively. Then the following holds true as τ tends to zero max kεn k2 = O(τ),

0≤n≤N

where the error εn = u(x, tn ) − un . 4

Proof. Subtracting (8) from (7) yields n X

εn+1 − εn = cτα

ωk ∆εn+1−k + Rn+1 .

(9)

k=0

Multiplying with εn+1 and integrating on domain Ω on both sides of equation (9), we arrive at 1 1 kεn+1 k22 − kεn k22 − kεn+1 k22 ≤ kεn+1 k22 − (εn , εn+1 ) 2 2 n X ωk (∆εn+1−k , εn+1 ) + (Rn+1 , εn+1 ). = cτα

(10)

k=0

Replacing n with p and summing up over p from 0 to n on both sides of (10), we obtain kεn+1 k2 ≤ 2cτα

p n X X p=0 k=0

α

= 2cτ

p n X X

n  X    ωk ∆ε p+1−k , ε p+1 + 2 R p+1 , ε p+1

ωk

p=0 k=0 p n X X α

= −2cτ

"

ωk

p=0 k=0

+2

p=0

2 p+1−k

∂ ε ∂x2

Ω

""

∂ε Ω

p+1−k

∂x

2 p+1−k

∂ ε + ∂y2 p+1

!

ε p+1 dxdy + 2

∂ε dxdy + ∂x

"

n  X  R p+1 , ε p+1 p=0

∂ε Ω

p+1−k

∂y

∂ε p+1 dxdy ∂y

#

n  X  R p+1 , ε p+1 p=0

    p  p " X n X  n X ∂ε p+1−k  ∂ε p+1 X ∂ε p+1−k  ∂ε p+1      dxdy   ωk = −2cτ +    ωk ∂x  ∂x ∂y  ∂y  Ω α

p=0

p=0

k=0

k=0

n  X  +2 R p+1 , ε p+1 p=0

≤ 2 max kε p k2 0≤p≤N

n X

kR p+1 k2

p=0

≤ 2C 0 max kε p k2 0≤p≤N

n X

τ2

p=0

p

0

≤ 2C max kε k2 (Nτ)τ 0≤p≤N

= C max kε p k2 τ. 0≤p≤N

In the third step of above derivation we use the integration by parts and the fact that εn vanishes on ∂Ω. Taking maximum over n yields max kεn k2 ≤ Cτ.

0≤n≤N

Assume that the solution of (8) satisfy uk (x)|∂Ω = 0, k = 0, 1, · · · . Totally similar to the proof of Theorem 2, we can also establish the priori estimation of the solution of (8) as following   n X   1 2 n+1 2 0 2 2 p+  2 (11) ku k2 ≤ C kψk2 + kϕk2 + τ kF k2  , n ≥ 0, p=0

0

where C is a constant independent of the time step τ. So the small perturbation in the initial conditions and source term would not cause large error in the numerical solution. Then we have the following theorem of stability. 5

Theorem 3. The semi-discretization (8) is unconditionally stable. 3. The MLS meshless method and full discretization The MLS meshless method is one of the most widely used meshless methods for its good properties of accuracy, robustness and high order of continuity. This section is mainly devoted to the construction of the shape function of the MLS meshless method. Suppose that the values ui = u(xi ) of the function u(x) to be determined is given at M node points xi , i = 1, 2, · · · , M. In the neighborhood Ω x of computation point x, u(x) can be approximated locally as [24, 35] uh (x, x) =

m X

pi (x)ai (x) = pT (x)a(x),

(12)

i=1

where x = [x, y]T is the space coordinates of point in the neighborhood Ω x of computation point x, pT (x) = [p1 (x), p2 (x), · · · , pm (x)], with pi (x) denoting the base function and m the number of base functions, aT (x) = [a1 (x), a2 (x), · · · , am (x)] is the coefficient vector to be determined. Usually we use monomial as base function. In special case we also have other alternatives such as singular or trigonometric function. The monomial base function in two dimensional space is linear base function: pT (x) = [1, x, y], m = 3, quadratic base function: pT (x) = [1, x, y, x2 , xy, y2 ], m = 6. Similarly in three dimensional space the monomial base function is linear base function: pT (x) = [1, x, y, z], m = 4, quadratic base function: pT (x) = [1, x, y, z, x2 , xy, y2 , yz, z2 , xz], m = 10. In the MLS approximation, the coefficients ai (x) are determined to assure that uh (x, x) is the best approximation for u(x) in the sense of weighted least squares in the neighborhood Ωx of computation point x. When discretizing the domain Ω with M node points, we can define a weight function Wi (x) = W(x − xi ) at each of the node points xi , i = 1, 2, · · · , M. The weight function Wi (x) is compactly supported, i.e. Wi (x) takes nonzero values in an infinite domain Ωi near xi and vanishes outside Ωi . We call Ωi the support domain of the weight function Wi (x) or of the point xi . Suppose that the neighborhood Ωx of computation point x contains N node points. The weighted summation of error of the approximation function uh (x, x) at these node points x = xi is  m 2 N N X h i2 X X  J= Wi (x) uh (x, xi ) − u(xi ) = Wi (x)  pk (xi )ak (x) − ui  . (13) i=1

i=1

k=1

To minimize J, we set

  m N X  ∂J X Wi (x)  =2 pk (xi )ak (x) − ui  p j (xi ) = 0, j = 1, 2, · · · , m. ∂a j (x) i=1 k=1

Then from (14) we obtain  N   N  m X X  X    Wi (x)pk (xi )p j (xi ) ak (x) =  Wi (x)p j (xi ) ui , j = 1, 2, · · · , m, k=1

i=1

(14)

(15)

i=1

namely

A(x)a(x) = B(x)u,

(16)

where A(x) =

N X

Wi (x)p(xi )pT (xi ),

i=1

6

  B(x) = W1 (x)p(x1 ) W2 (x)p(x2 ) · · · WN (x)p(xN ) ,

u = (u1 , u2 , · · · , uN )T .

We refer to ui as the node parameter. Because the number of node points contained in a support domain is usually larger than m the approximation function uh (x, x) does not pass through the nodal values [24]. From (16) the coefficient vector a(x) can be written as a(x) = A−1 (x)B(x)u.

(17)

Inserting (17) into (12) yields uh (x, x) = pT (x)A−1 (x)B(x)u = Φ(x, x)u,

(18)

where the shape function Φ(x, x) is formulated as Φ(x, x) = pT (x)A−1 (x)B(x).

(19)

The function u(x, x) is the locally best approximation of u(x) in the neighborhood Ωx of point x in the sense of weighted least squares. When the computation point x varies in the domain Ω, the evaluations of locally approximation function uh (x, x) at x = x form the global approximation function uh (x) of the function u(x) to be determined on the domain Ω, i.e. u(x) ≈ uh (x) = uh (x, x) x=x = Φ(x)u, (20)

where Φ(x) is

Φ(x) = Φ(x, x) x=x = pT (x)A−1 (x)B(x).

(21)

r = A−1 p.

(22)

In order to simplify the computation of derivatives of Φ(x), we let

Then r,i = A−1 (p,i − A,i r),

(23)

r,i j = A−1 (p,i j − A,i r, j − A, j r,i − A,i j r),

where the subscript “, i” represents the derivative with respect to x or y, while “, i j” represents the mixture derivative with respect to x and y. The first and second derivative of the shape function Φ(x) can be derived as following Φ,i = rT,i B + rT B,i ,

(24)

Φ,i j = rT,i j B + rT,i B, j + rT, j B,i + rT B,i j .

So far, there are still two significant issues left to be addressed. The first is determination of the dimension of support domain. Liu and Gu [24] suggest that generally circular or rectangular support domain can be used. The dimension (the radius when circular support domain is used, and the length and width when rectangular support domain is used) should be chosen to assure that the support domain contains sufficient number of node points to ensure the non-singularity of matrix A in (16). The support domain of a node point also should not be very large for the sake of cutting down the computation cost. To define the support domain for a node point xQ , the dimension of the support domain d s can be determined by [24] d s = α s dc , where α s is the dimensionless size of the support domain, and dc is the average nodal spacing near the point at xQ . Generally α s ∈ [2.0, 4.0] can yield good results. See [24] for the detail of obtaining dc . The second is the choose of 7

weight function. In practice, exponential function and spline function can be used as weight function. Liu and Gu [24] reported that the spline weight function is easier to use and handle than the exponential weight function. The most commonly used spline weight functions are the cubic spline weight function    ri ≤ 0.5, 2/3 − 4ri2 + 4ri3 ,     2 3 Wi (x) =  (25) 4/3 − 4ri + 4ri − 4ri /3, 0.5 < ri ≤ 1,     0, ri > 1, and the quartic spline weight function    1 − 6ri2 + 8ri3 − 3ri4 , 0 < ri ≤ 1, Wi (x) =   0, ri > 1,

(26)

where ri is the dimensionless distance from the computation point x to point xi , i.e. ri = kx − xi k/di , with di denoting the dimension of the support domain centered by xi . In the remainder of this section we exert ourselves to derive the full discretization for (8). Suppose that the problem domain is discretized by Nd internal node points and Nb boundary node points. We also assume that {xl |l = i1 , i2 , · · · , ini } denote the set of node points contained in the support domain of xi to construct the MLS shape function and Di = {i1 , i2 , · · · , ini }. Then we have the following Nd + Nb equations − cτα ω0 ∆uk+1 = uki + cτα uk+1 i i

k+1 X

k+1− j

ω j ∆ui

+ τψi + τFik+1 , i = 1, 2, · · · , Nd ,

(27)

j=1

uk+1 = g(xi , tk+1 ), i = Nd + 1, Nd + 2, · · · , Nd + Nb , i

(28)

where uij =

X

λlj Φ(i) l (xi ),

(29)

l∈Di

∆uij =

X

λlj ∆Φ(i) l (xi ).

(30)

l∈Di

The partial derivatives of shape function Φ can be obtained using (24). Let operator H = I − cτα ω0 ∆, with I representing the identity operator. Then we get the following full discretization for equation (8) at node point xi in each time step X

i λk+1 j HΦ j (xi ) =

j∈Di

X j∈Di

λkj Φij (xi ) + cτα

k+1 X l=1

ωl

X

λk+1−l ∆Φij (xi ) j

j∈Di

(31)

+ τψi + τFik+1 , i = 1, 2, · · · , Nd , X

i λk+1 j Φ j (xi ) = g(xi , tk+1 ), i = Nd + 1, Nd + 2, · · · , Nd + Nb .

(32)

j∈Di

It should be mentioned that equation (31) and (32) are system of nodal parameter λ and not of u itself. So we need to reconstruct u using (29) as λ is available. 4. Numerical examples To test the feasibility and efficiency of our approach and verify the theoretical results, two-dimensional numerical examples are solved with regular and irregular node points distribution in rectangular and circular domain. 8

Consider the following two-dimensional time-dependent fractional diffusion-wave equation in rectangular spatial domain Ω = [0, 1] × [0, 1] ∂α u(x, y, t) ∂2 u(x, y, t) ∂2 u(x, y, t) = + + f (x, y, t), (x, y) ∈ Ω, ∂tα ∂x2 ∂y2

(33)

u(x, y, 0) = 10x2 y2 (1 − x)(1 − y),

(34)

ut (x, y, 0) = 30x2 y2 (1 − x)(1 − y), (x, y) ∈ Ω, u(x, y, t) = 0, (x, y) ∈ ∂Ω,

(35)

where # t2−α t3−α + f (x, y, t) = 60x y (1 − x)(1 − y) Γ(4 − α) Γ(3 − α) i h − 20 (1 − 3x)y2 (1 − y) + (1 − 3y)x2 (1 − x) (t + 1)3 . 2 2

"

(36)

It can be easily checked that the exact solution of equation (33-36) is u(x, y, t) = 10x2 y2 (1 − x)(1 − y)(t + 1)3 .

(37)

In this section the error indicator we use is error =

M N 1 X X j ui − u(xi , t j ) , N × M i=1 j=1

(38) j

where N is as defined in section 2, M is number of spatial internal and boundary node points, and ui is the approximate solution obtained using the MLS meshless method. The convergence order is defined as co = log τ1

τ2

error(τ1 ) , error(τ2 )

where τ1 and τ2 are different time stepsizes, and error(τ) means that error corresponds to τ. In the following numerical experiments quartic base function and cubic spline weight function are used. The rectangular support domain is used and the dimensionless size α s of a support domain is set to be 3.0 for both the cases of rectangular and circular problem domain. 4.1. Numerical example with rectangular problem domain We first discretize the problem domain Ω = [0, 1] × [0, 1] with 3721 (61 × 61) regularly distributed points as shown in Fig.1 (a). Then the 2-D fractional diffustion-wave equation is solved using the presented meshless method and the results are compared with the exact solution. In Fig.2 the approximate and exact solution of equation (33-36) at t = 1 with the fractional order α = 1.8 are displayed. The absolute error is plotted in Fig.3. From Fig.2 and Fig.3 we see can see that the proposed meshless method is efficient and the results we obtain are reliable. We also discretize the problem domain Ω with 3721 irregularly distributed points (61 × 61 regularly distributed points then with small random perturbation in the coordinates) as shown in Fig.1 (b). The approximate and exact solution at t = 1 associated with irregularly distributed node points for α = 1.8 are plotted in Fig.4. The absolute error is also plotted in Fig.5. The approximation coincides well with the exact solution and the absolute error is almost of the same order of magnitude as the error for regular node points. Fig.6 is a curve of the average error versus time stepsize τ for α = 1.5, 1.8. In Table 1 we list the average error (defined in (38)) and convergence order of our approach for α = 1.5, 1.8, on rectangular domain with irregular node points distribution when the time stepsize τ decreases from 1/10 to 1/70. It can be seen that for both cases of α = 1.5 and 1.8 the average error decreases with the expected speed.

9

1 0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

y

y

1 0.9

0.4

0.4

0.3

0.3

0.2

0.2

0.1 0

0.1 0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

x

0.8

1

x

(a)

(b)

2

2

1.5

1.5 exact solution

approximate solution

Figure 1. Regular (a) and irregular (b) nodal distribution for rectangular domain

1

0.5

0 0

1

0.5 1 0.2

0.4

0 0

0.5 0.6

0.8

1

0

1 0.2

0.4

0.5 0.6

0.8

y

x

1

0

y

x

(a)

(b)

Figure 2. Approximate and exact solution at t = 1 in rectangular domain with α = 1.8, τ = 1/70, for regular node points distribution

0.015

e

0.01

0.005

0 0

1 0.2

0.4

0.5 0.6

0.8

1

0

y

x

Figure 3. Absolute value of the difference between approximate and exact solution at t = 1 with α = 1.8, τ = 1/70, for regular node points distribution

10

2

1.5

1.5 exact solution

approximate solution

2

1

0.5

0 0

1

0.5 1 0.2

0.4

1

0 0

0.2

0.5 0.6

0.8

1

0

0.5

0.4

0.6

0.8

y

x

1

0

y

x

(a)

(b)

Figure 4. Approximate and exact solution at t = 1 in rectangular domain with α = 1.8, τ = 1/70, for irregular node points distribution

0.015

e

0.01

0.005

0 0

1 0.2

0.4

0.5 0.6

0.8

1

0

y

x

Figure 5. Absolute value of the difference between approximate and exact solution at t = 1 with α = 1.8, τ = 1/70, for irregular node points distribution

11

0.022 α=1.5 α=1.8

0.02 0.018 0.016

error

0.014 0.012 0.01 0.008 0.006 0.004 0.002

0.1

0.08

0.06 τ

0.04

0.02

Figure 6. Error as a function of time stepsize τ at t = 1, with α = 1.5, 1.8, for irregular node points distribution shown in Fig.1 (b)

Table 1 Average error and convergence order associated with different time stepsize for irregular node points distribution on rectangular domain (α = 1.5, 1.8, t = 1)

α = 1.5 stepsize τ 1/10 1/20 1/30 1/40 1/50 1/60 1/70

error 0.0125502487 0.0066277363 0.0045292253 0.0034518685 0.0027951743 0.0023526010 0.0020338938

α = 1.8 stepsize τ 1/10 1/20 1/30 1/40 1/50 1/60 1/70

convergence order 0.9211 0.9390 0.9442 0.9457 0.9454 0.9443

error 0.0204967111 0.0106960039 0.0072811829 0.0055407991 0.0044822148 0.0037688513 0.0032548236

convergence order 0.9383 0.9485 0.9495 0.9502 0.9508 0.9512

4.2. Numerical example with circular problem domain In order to test the ability of the proposed method to handle problem on complex-shaped domain, we solve the 2-D FDWE with the MLS meshless approach on circular problem domain. Specifically the equation we solve is equation (33-36) with (35) slightly modified to n assure that the boundary conditiono is just the value of exact solution on the boundary. First the circular domain Ω = (x, y) (x − 0.5)2 + (y − 0.5)2 ≤ 0.25 is discretized with 3721 irregularly distributed node points as shown in Fig.7. Then we solve equation (33-36) with the MLS meshless approach using these irregularly distributed node points for both α = 1.5 and 1.8.

12

1 0.9 0.8 0.7

y

0.6 0.5 0.4 0.3 0.2 0.1 0

0.2

0

1

0.8

0.6

0.4 x

Figure 7. Irregular nodal distribution for circular domain

2

2

1.5

1.5

exact solution

approximate solution

Fig.8 is a comparison of the approximate solution with exact solution on Ω for α = 1.5. In Fig.9 the absolute error of our approach at t = 1 with α = 1.5 is also plotted. These two figures indicate the reliability of the outputs of our code and the efficiency of the proposed MLS meshless method. Fig. 10 is the curve of the average error as a function of time stepsize τ when t = 1, α = 1.5, 1.8. Similar to the previous subsection, we list the average error and convergence order for α = 1.5, 1.8, as shown in Table 2 from which we can see that the proposed method also works well on circular domain.

1

0.5

0 0

1

0.5 1 0.2

0.4

0 0

0.5 0.6

0.8

1

0

1 0.2

0.4

0.5 0.6

0.8

y

x

1

0

y

x

(a)

(b)

Figure 8. Approximate and exact solution at t = 1 in circular domain with α = 1.5, τ = 1/70, for irregular node points distribution

13

−3

x 10 8

e

6

4

2

0 0

1 0.2

0.4

0.5 0.6

0.8

1

0

y

x

Figure 9. Absolute value of the difference between approximate and exact solution at t=1 with α = 1.5, τ = 1/70, for irregular node points distribution on circular domain

0.02 α=1.5 α=1.8

0.018 0.016 0.014

error

0.012 0.01 0.008 0.006 0.004 0.002 0

0.1

0.08

0.06 τ

0.04

0.02

Figure 10. Error as a function of time stepsize τ at t=1.0, with α = 1.5, 1.8, for irregular node points distribution shown in Fig.7

Table 2 Average error and convergence order associated with different time stepsize for irregular node points distribution on circular domain (α = 1.5, 1.8, t = 1)

α = 1.5 stepsize τ 1/10 1/20 1/30 1/40 1/50 1/60 1/70

error 0.0116516710 0.0060965285 0.0041315979 0.0031234079 0.0025090268 0.0020949923 0.0017968573

α = 1.8 stepsize τ 1/10 1/20 1/30 1/40 1/50 1/60 1/70

convergence order 0.9345 0.9595 0.9724 0.9816 0.9892 0.9958

error 0.0189309749 0.0097300702 0.0065580462 0.0049523355 0.0039794597 0.0033253801 0.0028547748

convergence order 0.9602 0.9730 0.9762 0.9801 0.9849 0.9899

5. Conclusions In this paper, we first transform the original 2-D fractional diffusion-wave equation into an equivalent integrodifferential equation and then semi-discretize the transformed equation. We approximate the integration term with 14

Lubich’s first order fractional multistep method and discuss the convergence and stability properties of the semidiscretization. Then we obtain a full discretization for 2-D FDWE using the MLS meshless approach. Finally we verify the theoretical results by numerical examples on rectangular and circular domain with regular and irregular node points distribution. In addition we find that large α leads to relatively large average error. The reason is that large value of α probably leads to large error constant in the truncation error. In summary, the presented meshless method is an efficient and promising approach for 2-D FDWE. It makes the problem on domain with complex shape easier to be handled than FDM and has a good potential of use in engineering and science. Acknowledgements The authors owe thanks to the financial support of National Natural Science Foundation of China (Grant Nos. 60931002, 11001072 and 11101381). References References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional integrals and derivatives: theory and applications, Gordon and Breach, Yverdon, 1993. K. Diethelm, The analysis of fractional differential equations, Springer, Berlin, 2010. I. Podlubny, Fractional differential equations, Academic Press, New York, 1999. T.S. Aleroev, H.T. Aleroeva, J.F. Huang, N.M. Nie, Y.F. Tang, S.Y. Zhang, Features of seepage of a liquid to a chink in the cracked deformable layer, Int. J. Modell. Simul. Sci. Comput. 1 (3) (2010) 333-347. R.R. Nigmatullin, To the theoretical explanation of the universal response, Physica Status (B): Basic Res. 123 (2) (1984) 739-745. R.R. Nigmatullin, Realization of the generalized transfer equation in a medium with fractal geometry, Physica Status (B): Basic Res. 133 (1) (1986) 425-430. W. Wess, The fractional diffusion equation, J. Math. Phys. 27 (11) (1996) 2782-2785. W.R. Schneider, W. Wess, Frational diffusion and wave equations, J. Math. Phys. 30 (1) (1989) 134-144. F. Mainardi, Frational relaxation-oscilation and fractional diffusion-wave phenomena, Chaos Solitons Fractals 7 (9) (1996) 1461-1477. F. Mainardi, The fundamental solutions for the fractional diffusion-wave equation, Appl. Math. Lett. 9 (6) (1996) 23-28. O.P. Agrawal, A general solution for the fourth-order fractional diffusion-wave equation, Fract. Calculus. Appl. Anal. 3 (1)(2000) 1-12. O.P. Agrawal, A general solution for the fourth-order fractional diffusion-wave equation defined in bounded domain. Comput. Struct. 79 (2001) 1497-1501. O.P. Agrawal, Response of a diffusion-wave system subjected to deterministic and stochastic fields, Z. Angew. Math. Mech. 83 (4) (2003) 265-274. R. Metzler, J. Klafler, Boundary value problems for fractional diffusion equation, Physica A 278 (2000) 107-125. S.B. Yuste, L. Acedo, An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equation, SIAM J. Numer. Anal. 42 (5) (2005) 1862-1874. S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (1) (2006) 264-274. F. Liu, S. Chen, V. Anh, I. Turner, Analysis of a discrete non-Markovian random walk approximation for the time fractional diffusion equation, ANZIAM J. 46 (E) (2005) 488-504. P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2) (2008) 1079-1095. X.J. Li, C.J. Xu, A space-time spectral method for the time fractional differential equation, SIAM J. Numer. Anal. 47 (3) (2009) 2108-2131. Y.M. Lin, C.J. Xu, Finite difference/spectral approximation for the time-fractional diffusion equation, J. Comput. Phys. 225 (2007) 1533-1552. Z.Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209. C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205-213. R.A. Gingold, J.J. Moragham, Smooth particle hydrodynamics: theory and applications to non spherical stars, Monthly Notices of the Royal Astronomical Society 181 (1977) 375-389. G.R. Liu, Y.T. Gu, An introduction to meshfree methods and their programming, Springer, Berlin, 2005. E.J. Kansa, Multiquadrics-A scattered data approximation scheme with applications to computational fluid dynamics, Comput. Math. Appl. 19 (8/9) (1990) 127-145. E. Onate, S. Idelsohn, O.C. Zienkiewicz, R.L. Taylor, C. Sacco, A finite method in computational mechanics: applications to convective transport and fluid flow, Int. J. Numer. Meth. Eng. 39 (1996) 3839-3866. Y.T. Gu, Meshfree methods and their comparisons, Int. J. Comput. Meth. 2(4) (2005) 477-515. Y.T. Gu, P. Zhuang, F. Liu, An advanced implicit meshless approach for the non-linear anomalous subdiffusion equation, Comput. Modell. Eng. Sci. 56 (3) (2010) 303-333. P. Zhuang, Y.T. Gu, F. Liu, I. Turner, P.K.D.V. Yarlagadda, Time-dependent fractional advection-diffusion equations by an implicit MLS meshless method, Int. J. Numer. Meth. Eng. 88 (2011) 1346-1362.

15

[30] Q. Liu, Y.T. Gu, P. Zhuang, F. Liu, Y. Nie, An implicit RBF meshless approach for time fractional diffusion equations, Comput. Mech. 48 (2011) 1-12. [31] Q. Liu, F. Liu, I. Turner, V. Anh, Y.T. Gu, A RBF mesheless approach for modeling a fractal mobile/immobile transport model, Appl. Math. Comput. 226 (2014) 336-347. [32] J.F. Huang, Y.F. Tang, V. Luis, J.Y. Yang, Two finite difference schemes for time fractional diffusion-wave equation. Numer. Algor. 64 (4) (2013) 707-720. [33] Ch. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (3) (1986) 704-719. [34] J.C. Lopze-Marcos, A difference scheme for a nonlinear partial integrodifferential equation, SIAM J. Numer. Anal. 27 (1) (1990) 20-31. [35] X. Zhang, Y. Liu, Meshless methods (in Chinese), Tsinghua University Press, Beijing, 2004.

16