Finite Difference and Sinc-Collocation Approximations to a Class of Fractional Diffusion-Wave Equations

We propose an efficient numerical method for a class of fractional diffusion-wave equations with the Caputo fractional derivative of order α. This approach is based on the finite difference in time and the global sinc collocation in space. By utilizing the collocation technique and some properties of the sinc functions, the problem is reduced to the solution of a system of linear algebraic equations at each time step. Stability and convergence of the proposedmethod are rigorously analyzed.Thenumerical solution is of 3 − αorder accuracy in time and exponential rate of convergence in space. Numerical experiments demonstrate the validity of the obtained method and support the obtained theoretical results.


Introduction
Fractional diffusion-wave equations have attracted considerable attention as generalization of the classical diffusion/wave equation by replacing the integer-order time derivative with a fractional derivative of order  (1 <  < 2) [1].It can be derived from the physical case of a general random timedependent velocity function equipped with an algebraic timecorrelation function (but with a Gaussian space correlation) [2].Owing to the anomalous superdiffusion of particles, many electromagnetic, acoustic, mechanical, and biological responses can be modelled accurately by fractional diffusionwave equations, for example, the propagation of stress waves in viscoelastic solids [3,4].There exist some analytical methods to find exact solutions of fractional diffusion-wave equations [5][6][7][8][9].However, it is usually difficult or even impossible to achieve the exact solutions in many cases, and ones have to resort to numerical methods.
Compared with the considerable literature on numerical solutions of fractional diffusion equations [10][11][12][13][14][15], only a few works have been carried out in the numerical methods of fractional diffusion-wave equations.The homotopy analysis method [16] and the Adomian decomposition method [17] were used to construct analytical approximate solutions of fractional diffusion-wave equations, respectively.Sun and Wu [18] presented a full discrete difference scheme for the initial-boundary-value problem of a diffusion-wave equation by introducing two new variables and obtained the stability and convergence properties by using energy method.In [19], under the weak smoothness conditions, two finite difference schemes with first-order accuracy in temporal direction and second-order accuracy in spatial direction were proposed to solve the same problem.Murillo and Yuste [20] developed an explicit difference method where the L2 discretization formula is used for solving fractional diffusionwave equations.Due to the nonlocal nature of fractional derivatives, all previous solutions have to be saved to compute the solution at the current time level.It makes the storage very expensive when low-order methods are employed for spatial discretization.
In this paper, we propose a numerical method mixing finite difference with sinc collocation to solve a class of fractional diffusion-wave equations.The sinc method is widely used to solve integral, integro-differential, ordinary differential, and partial differential equations [21][22][23][24][25][26][27].As far as we know, there are very limited papers on solving fractional differential equation by the sinc method [28,29].The present method has great advantage in storage requirement because the sinc-collocation method needs fewer grid points to produce highly accurate solution when it is compared with a low-order method.
The remainder of this paper is organized as follows.In Section 2, we introduce some necessary definitions and relevant results for developing this method.Section 3 is devoted to the finite difference and sinc-collocation discretization for the fractional diffusion-wave equations.As a result, a system of linear algebraic equations is formed, and the solutions of the considered problems are obtained.Section 4 is concerned with the stability and convergence analysis.In Section 5, some numerical experiments are given to demonstrate the effectiveness of the proposed method and the obtained theoretical results.A brief conclusion ends this paper in the final section.

Notations and Some Preliminary Results
In this section, we introduce some basic definitions and relevant results of the fractional calculus [30,31] and sinc functions.
Definition 2 (see [32]).Let  ∈ R + and  = ⌈⌉.The Caputo fractional differential operator     for  ≤  ≤  is defined as The sinc functions and its properties are discussed thoroughly in [21,23].For any ℎ > 0, the translated sinc functions with equidistant space nodes are given as where the sinc functions are defined on the whole real line by If  is defined on R, then for any ℎ > 0 the series is called the Whittaker cardinal expansion of  whenever this series converges. can be approximated by truncating (5).To construct our needed approximations on the interval [, ], we choose which maps a finite interval [, ] onto R.  is a conformal map which maps the eye-shaped domain in the complex plane onto the infinite strip   of the complex z-plane The basis functions on [, ] are taken to be the composite translated sinc functions Thus one may define the inverse image of the equidistant space node {ℎ} as The class of functions such that the known exponential convergence rate exists for the sinc interpolation is denoted by (  ) and defined in the following.
The following theorem gives the error which results from differentiating the truncated cardinal series.
Theorem 4 (see [21,22] for all  = 0, 1, . . ., ; here  3 is a constant depending only on , , , , , and . The above expressions show that the sinc interpolation on (  ) converges exponentially.We also require the following derivatives of the composite translated sinc functions evaluated at the nodes:

The Finite Difference and Sinc-Collocation Discretization
We consider the following fractional diffusion-wave equations with a nonhomogeneous field [33]: with the initial conditions and the boundary conditions where  ∈ [, ] and  > 0 are space and time variables,  and  are constants, and (, ) denotes the field variable.
Here the time-fractional derivative is defined as the Caputo fractional derivative, and 1 <  < 2. Many authors refer to the fractional equation ( 19) as the fractional diffusion-wave equation when 1 <  < 2, which is expected to interpolate the diffusion equation and the wave equation [7,18].

Temporal Discretization by a Finite Difference Scheme.
First, we derive a finite difference scheme for temporal discretization of this equation.Let   := ,  = 0, 1, 2, . .., where  is the time stepsize.In order to discretize the time-fractional derivative by using a finite difference approximation [18,34], we introduce the following lemmas.
A rigorous analysis of the convergence rate will be provided later.The above scheme can be rewritten into Specially, for the case  = 1, the scheme simply reads So ( 37) and (38) together with the initial condition and the boundary conditions form a complete semidiscrete problem.

Space Discretization by the Sinc-Collocation Method.
Next we consider space discretization to (37) by the sinccollocation method.We select the collocation points   by (10).The space discretization proceeds by approximating the solution based on the composite translated sinc functions (9) The unknown parameters    will be determined by the collocation method.It is noted that the approximation in (41) satisfies the boundary conditions in (40) since   (, ℎ)(),  = −, − + 1, . . .,  are zero when  tends to  and .Now substituting (41) into (37) and collocating in 2 + 1 points   , we obtain where  = −, − + 1, . . ., .Based on ( 17) and ( 18), we let Hence, (42) is reduced immediately by ( 16) and (43) to where  = −, − + 1, . . ., .To obtain a matrix representation of the above equation, we let where the matrix  is the identity matrix and  = 2, 3, . . ., −1.Therefore, at each time step, we get the following system of 2 + 1 linear equations with 2 + 1 unknown parameters    , and this system can be expressed in a matrix form where For the initial condition (39), we have Consequently (46) can be solved easily for the unknown coefficients U  .Hence the approximation solution given in (41) can be obtained.

Stability and Convergence Analysis of the Derived Method
To analyse the stability and convergence of the derived method, the inner product is defined by ⟨, ⟩ = ∫     with the corresponding norm ‖‖ = (, ) (1/2) , which will be used thereafter.We first give the following lemmas.
(54) Applying Lemma 6 to the term in left side of (54), we obtain On the other side, based on the fact that
where  is constant.
Remark 10.In Theorem 9, the error estimate formula contains  2  −  (3/2)  −√ (the second term in the right-hand side), where the error contribution from the spatial approximation is affected by the inverse of the time step.It is worthwhile noting that similar results are also found for the classical diffusion/wave equation.However for large ,  (3/2)  −√ can be much smaller than ; therefore, this affection generally would not reduce the global accuracy.

Numerical Examples
To validate the effectiveness of the proposed method for the problem ( 19)-( 21), we consider the example given in [18].Consider the following: The exact solution of the above problem is [33] where   () = ∑ ∞ =0   /Γ( + 1) which is one-parameter Mittag-Leffler function.
For solving the above problem (with  = 1.3 and 1.7, resp.) by using the method described in Section 3, we choose  = 1 and  = /2, and this leads to ℎ = / √ 2.We will report the efficiency and accuracy of the given method based on the  2 -errors and  ∞ -errors.Figure 1   diagrams of the numerical solutions   (, ) on the whole computational domain [0, 1] × [0, 1] with  = 0.001,  = 20.In Figure 2, we plot the curves of the numerical solutions   (, ) and exact solutions (, ) for several fixed time instants.A good agreement of the numerical solution with the exact one is achieved.Furthermore, Figure 3 shows the absolute error |  (, ) − (, )| obtained by the present method with  = 0.001,  = 20.
To check the convergence behavior of the numerical solution with respect to the time step  and the parameter  about the number of collocation points, we represent the errors ‖  −    ‖ in two discrete norms:  ∞ and  2 .All the numerical results reported in Figures 4 and 5, and Tables 1 and 2 are evaluated at  = 0.8.
Firstly, the computational investigation is concerned with the temporal convergence rate.We choose the parameter  = 60 and it is a value large enough such that the errors coming from the spatial approximation are negligible [35].In Tables 1 ( = 1.3) and 2 ( = 1.7), we list the temporal errors when  decreases from 1/10 to 1/320 and the convergence order which is very close to the expectation order 3 − .We also plot the errors in the  ∞ and  2 norms as a function of the time stepsize  for several  in Figure 4, where a logarithmic scale has been used for both -axis and error axis.From Figure 4, it is clear that, for  = 1.3 and 1.7, the slopes of the error curves in these log-log plots approach 1.7 and 1.3, respectively.So the proposed method yields a temporal approximation order which is close to 3− as forecasted by the theoretical estimate.

Figure 5 :
Figure 5: Spatial errors in the  ∞ and  2 norms versus  for several .