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

Abstract

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 proposed method are rigorously analyzed. The numerical 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.

1. 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 time-dependent velocity function equipped with an algebraic time-correlation 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 diffusion-wave 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–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–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 diffusion-wave 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–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.

2. 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 1 (see [32]). Let 𝛼∈𝓡₊. The operator 𝐽^𝛼_𝑎 defined on 𝐿₁[𝑎,𝑏] by 𝐽^𝛼_𝑎𝑓(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝑠)^𝛼−1𝑓(𝑠)𝑑𝑠 (1) for 𝑎≤𝑡≤𝑏 is called the Riemann-Liouville fractional integral operator of order 𝛼. For 𝛼=0, one sets 𝐽⁰_𝑎∶=𝐼, that is, the identity operator.

Definition 2 (see [32]). Let 𝛼∈𝓡₊ and 𝑛=⌈𝛼⌉. The Caputo fractional differential operator ^𝐶𝐷^𝛼_𝑎 for 𝑎≤𝑡≤𝑏 is defined as ^𝐶𝐷^𝛼_𝑎𝑓(𝑡)=𝐽_𝑎^𝑛−𝛼𝐷^𝑛𝑓(𝑡)=1Γ(𝑛−𝛼)^𝑡_𝑎(𝑡−𝑠)^{𝑛−𝛼−1}𝑓^(𝑛)(𝑠)𝑑𝑠. (2)

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 𝑆(𝑘,ℎ)(𝑧)=sinc𝑧−𝑘ℎℎ,𝑘=0,±1,±2,…, (3) where the sinc functions are defined on the whole real line by (𝑥)=()sincsin𝜋𝑥𝜋𝑥,𝑥≠0,1,𝑥=0. (4) If 𝑓 is defined on ℝ, then for any ℎ>0 the series 𝐶𝑓,ℎ(𝑧)=^∞_𝑘=−∞𝑓(𝑘ℎ)𝑆(𝑘,ℎ)(𝑧) (5) 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 𝜙(𝑥)=ln𝑥−𝑎𝑏−𝑥 (6) which maps a finite interval [𝑎,𝑏] onto ℝ. 𝜙 is a conformal map which maps the eye-shaped domain in the complex 𝑧-plane 𝐷_𝐸=𝑧∈ℂ∶arg𝑧−𝑎𝜋𝑏−𝑧<𝑑≤2 (7) onto the infinite strip 𝐷_𝑆 of the complex 𝑧-plane 𝐷_𝑆=|(𝑧𝜋𝑧∈ℂ∶Im)|<𝑑≤2. (8)

The basis functions on [𝑎,𝑏] are taken to be the composite translated sinc functions 𝑆_𝜙(𝑥)()𝜙(𝑥)𝜙(𝑥)𝑘,ℎ)(=𝑆𝑘,ℎ=sinc−𝑘ℎℎ. (9) Thus one may define the inverse image of the equidistant space node {𝑖ℎ} as 𝑥_𝑖=𝜙⁻¹(𝑖ℎ)=𝑎+𝑏𝑒^𝑖ℎ1+𝑒^𝑖ℎ,𝑖=0,±1,±2,…. (10) The class of functions such that the known exponential convergence rate exists for the sinc interpolation is denoted by 𝐵(𝐷_𝐸) and defined in the following.

Definition 3 (see [21, 23]). Let 𝐵(𝐷_𝐸) be the class of functions 𝑓 which are analytic in 𝐷_𝐸 and satisfy _{𝜙−1(𝑥+𝐿)}𝑓(𝑧)𝑑𝑧⟶0,𝑥⟶±∞, (11) where 𝐿={𝑖𝜐∶|𝜐|<𝑑≤(𝜋/2)} and _𝜕𝐷𝐸𝑓(𝑧)𝑑𝑧<∞ (12) on the boundary of 𝐷_𝐸 (denoted 𝜕𝐷_𝐸).

The following theorem gives the error which results from differentiating the truncated cardinal series.

Theorem 4 (see [21, 22]). Let 𝜙^′𝑓/𝑔∈𝐵(𝐷_𝐸) and let sup_{𝑥∈[𝑎,𝑏],−𝜋/ℎ≤𝑡≤𝜋/ℎ}dd𝑥^𝑛𝑔(𝑥)𝑒^{𝑖𝑡𝜙(𝑥)}≤𝐶₁ℎ^−𝑛, (13) for 𝑛=0,1,…,𝑚; here 𝑔 is a weight function and 𝐶₁ is a constant depending only on 𝑚, 𝜙, and 𝑔. Assume that there exists a constant 𝐶₂ such that 𝑓(𝑥)𝑔(𝑥)≤𝐶₂𝑒^{𝛽𝜙(𝑥)}1+𝑒^{𝜙(𝑥)2𝛽}; (14) here 𝛽 is a positive constant. Then taking ℎ=𝜋𝑑/𝛽𝑁, one has d^𝑛𝑓(𝑥)d𝑥^𝑛−^𝑁_𝑘=−𝑁𝑓𝑥_𝑘𝑔𝑥_𝑘d^𝑛d𝑥^𝑛𝑔(𝑥)𝑆_𝜙(𝑘,ℎ)(𝑥)≤𝐶₃𝑁^(𝑚+1)/2exp−𝜋𝑑𝛽𝑁 (15) for all 𝑛=0,1,…,𝑚; here 𝐶₃ 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: 𝛿⁽⁰⁾_𝑘,𝑖=𝑆_𝜙(𝑘,ℎ)(𝑥)_{𝑥=𝑥𝑖}=1,𝑘=𝑖,0,𝑘≠𝑖.𝛿⁽¹⁾_𝑘,𝑖=d𝑆d𝜙^𝜙(𝑘,ℎ)(𝑥)_{𝑥=𝑥𝑖}=()0,𝑘=𝑖,−1^𝑖−𝑘(𝑖−𝑘)ℎ,𝑘≠𝑖.𝛿⁽²⁾_𝑘,𝑖=d²d𝜙²𝑆_𝜙(𝑘,ℎ)(𝑥)_{𝑥=𝑥𝑖}=−𝜋²3ℎ²(),𝑘=𝑖,−2−1^𝑖−𝑘(𝑖−𝑘)²ℎ²,𝑘≠𝑖. (16)(17)(18)

3. The Finite Difference and Sinc-Collocation Discretization

We consider the following fractional diffusion-wave equations with a nonhomogeneous field [33]: 1𝑐𝜕^𝛼𝑢(𝑥,𝑡)𝜕𝑡^𝛼=𝜕²𝑢(𝑥,𝑡)𝜕𝑥²+1𝐾𝑓(𝑥,𝑡),𝑎<𝑥<𝑏,𝑡>0 (19) with the initial conditions 𝑢()(𝑥),()𝑥,0=𝜑𝜕𝑢𝑥,0(𝑥)𝜕𝑡=𝜓,𝑎<𝑥<𝑏 (20) and the boundary conditions 𝑢(𝑎,𝑡)=0,𝑢(𝑏,𝑡)=0,𝑡>0, (21) 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].

3.1. 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.

Lemma 5 (see [18]). Suppose 𝑓(𝑡)∈𝐶²[0,𝑡_𝑛]. Then ^𝑡𝑛₀𝑓^′(𝑡)𝑑𝑡𝑡_𝑛−𝑡^𝛼−1−𝜏^1−𝛼×𝑏2−𝛼₀𝑓𝑡_𝑛−^𝑛−1_𝑘=1𝑏_𝑘−1−𝑏_𝑘𝑓𝑡_𝑛−𝑘−𝑏_𝑛−1𝑓𝑡₀≤12−𝛼2−𝛼+212^3−𝛼−3−𝛼1+2^1−𝛼max_{0≤𝑡≤𝑡𝑛}𝑓^′′(𝑡)𝜏^3−𝛼, (22) where 1<𝛼<2 and 𝑏_𝑘=(𝑘+1)^2−𝛼−𝑘^2−𝛼,𝑘=0,1,2,…. (23)

It is direct to check that 1=𝑏₀>𝑏₁>⋯>𝑏_𝑛>0,𝑏_𝑛⟶0as𝑛⟶∞. (24)

Let 𝜙()=𝜕𝑥,𝑡^𝛼𝑢(𝑥,𝑡)𝜕𝑡^𝛼=1Γ(2−𝛼)^𝑡₀𝜕²𝑢(𝑥,𝑠)𝜕𝑠²𝑑𝑠(𝑡−𝑠)^𝛼−1, (25) and define the temporal grid functions 𝑈^𝑘=𝑢𝑥,𝑡_𝑘,Φ^𝑘=𝜙𝑥,𝑡_𝑘,𝑘=0,1,2,…. (26) By the Taylor expansion, it follows that 𝑈_𝑡^𝑘−(1/2)=1𝜏𝑈^𝑘−𝑈^𝑘−1+𝑟^𝑘₁𝜏²,1𝑐Φ^𝑘−(1/2)=𝑈^𝑘−(1/2)_𝑥𝑥+1𝐾𝑓^𝑘−(1/2)+𝑟^𝑘₂𝜏², (27)(28) where 𝑟^𝑘₁, 𝑟^𝑘₂ are constants. Based on Lemma 5, one has Φ^𝑘=𝜏^1−𝛼Γ()𝑏3−𝛼⁰𝑈^𝑘_𝑡−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑈_𝑡^𝑘−𝑙−𝑏_𝑘−1𝑈⁰_𝑡+𝑂𝜏^3−𝛼,𝑘=1,2,3,…. (29) Hence, Φ^𝑘−(1/2)=12Φ^𝑘+Φ^𝑘−1=𝜏^1−𝛼Γ()×𝑏3−𝛼₀𝑈_𝑡^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑈_𝑡^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝑈⁰_𝑡+𝑟^𝑘₃𝜏^3−𝛼, (30) where 𝑟^𝑘₃ are constants. Let 𝛿_𝑡𝑈^𝑘−(1/2)=(1/𝜏)(𝑈^𝑘−𝑈^𝑘−1). In view of (27), one has Φ^𝑘−(1/2)=𝜏^1−𝛼Γ()×𝑏3−𝛼₀𝛿_𝑡𝑈^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑈^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝑈⁰_𝑡+𝜏^3−𝛼Γ()𝑏3−𝛼⁰𝑟^𝑘₁−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑟₁^𝑘−𝑙+𝑟^𝑘₃𝜏^3−𝛼. (31) Then substituting (31) into (28) and observing 𝑈⁰_𝑡=𝜓(𝑥), one obtains 𝜏^1−𝛼()𝑏𝑐Γ3−𝛼⁰𝛿_𝑡𝑈^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑈^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜓=𝑈^𝑘−(1/2)_𝑥𝑥+1𝐾𝑓^𝑘−(1/2)+𝑅^𝑘_𝑇, (32) where 𝑅^𝑘_𝑇𝜏=−^3−𝛼()𝑏𝑐Γ3−𝛼⁰𝑟^𝑘₁−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑟₁^𝑘−𝑙−1𝑐𝑟^𝑘₃𝜏^3−𝛼+𝑟^𝑘₂𝜏². (33)

It follows from (32) that one can construct easily the following semidiscrete finite difference scheme for (19): 𝜏^1−𝛼()𝑏𝑐Γ3−𝛼⁰𝛿_𝑡𝑢^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑢^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜓=𝑢^𝑘−(1/2)_𝑥𝑥+1𝐾𝑓^𝑘−(1/2), (34) which is equivalent to 1𝑐𝜏^𝛼Γ()Δ3−𝛼¹+Δ₂=12𝑢^𝑘_𝑥𝑥+𝑢^𝑘−1_𝑥𝑥+1𝐾𝑓^𝑘−(1/2),𝑘=2,3,…, (35) where Δ₁=𝑏₀𝑢^𝑘−2𝑏₀−𝑏₁𝑢^𝑘−1,Δ₂=^𝑘−1_𝑙=2𝑏_𝑙−2−2𝑏_𝑙−1+𝑏_𝑙𝑢^𝑘−𝑙+𝑏_𝑘−2−𝑏_𝑘−1𝑢⁰−𝑏_𝑘−1𝜏𝜓. (36) According to (33), this scheme is of (3−𝛼)-order accuracy. A rigorous analysis of the convergence rate will be provided later. The above scheme can be rewritten into 𝑏₀𝑐𝜏^𝛼Γ()𝑢3−𝛼^𝑘−12𝑢^𝑘_𝑥𝑥=2𝑏₀−𝑏₁𝑐𝜏^𝛼Γ()𝑢3−𝛼^𝑘−1−1𝑐𝜏^𝛼Γ()Δ3−𝛼²+12𝑢^𝑘−1_𝑥𝑥+1𝐾𝑓^𝑘−(1/2). (37) Specially, for the case 𝑘=1, the scheme simply reads 𝑏₀𝑐𝜏^𝛼Γ()𝑢3−𝛼¹−12𝑢¹_𝑥𝑥=𝑏₀𝑐𝜏^𝛼Γ()𝑢3−𝛼⁰+1+𝜏𝜓2𝑢⁰_𝑥𝑥+1𝐾𝑓^(1/2). (38) So (37) and (38) together with the initial condition 𝑢⁰=𝑢(𝑥,0)=𝜑(𝑥),𝑎<𝑥<𝑏 (39) and the boundary conditions 𝑢^𝑘_𝑎=𝑢𝑎,𝑡_𝑘=0,𝑢^𝑘_𝑏=𝑢𝑏,𝑡_𝑘=0,𝑘≥1 (40) form a complete semidiscrete problem.

3.2. Space Discretization by the Sinc-Collocation Method

Next we consider space discretization to (37) by the sinc-collocation method. We select the collocation points 𝑥_𝑖 by (10). The space discretization proceeds by approximating the solution based on the composite translated sinc functions (9) 𝑢^𝑘≃𝑢^𝑘_𝑁=^𝑛_𝑗=−𝑛𝑢^𝑘_𝑗𝑆_𝜙𝑗,ℎ(𝑥). (41) 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 𝑏₀𝑐𝜏^𝛼Γ(3−𝛼)^𝑛_𝑗=−𝑛𝑢^𝑘_𝑗𝑆_𝜙𝑗,ℎ𝑥_𝑖−12^𝑛_𝑗=−𝑛𝑢^𝑘_𝑗𝑆_𝜙^′′𝑗,ℎ𝑥_𝑖=2𝑏₀−𝑏₁𝑐𝜏^𝛼Γ(3−𝛼)^𝑛_𝑗=−𝑛𝑢_𝑗^𝑘−1𝑆_𝜙𝑗,ℎ𝑥_𝑖+12^𝑛_𝑗=−𝑛𝑢_𝑗^𝑘−1𝑆_𝜙^′′𝑗,ℎ𝑥_𝑖−1𝜏^𝛼Γ()×3−𝛼_{𝑘−1𝑙=2}𝑏_𝑙−2−2𝑏_𝑙−1+𝑏_𝑙×^𝑛_𝑗=−𝑛𝑢_𝑗^𝑘−𝑙𝑆_𝜙𝑗,ℎ𝑥_𝑖+𝑏_𝑘−2−𝑏_𝑘−1𝑢⁰𝑥_𝑖−𝜏𝑏_𝑘−1𝜓𝑥_𝑖+1𝐾𝑓^𝑘−(1/2)𝑥_𝑖, (42) where 𝑖=−𝑛,−𝑛+1,…,𝑛. Based on (17) and (18), we let 𝑞_𝑗𝑖∶=𝑆_𝜙^′′𝑗,ℎ𝑥_𝑖=𝜙^′′𝑥_𝑖𝛿⁽¹⁾_𝑗,𝑖+𝜙^′𝑥_𝑖²𝛿⁽²⁾_𝑗,𝑖. (43) Hence, (42) is reduced immediately by (16) and (43) to 𝑏₀𝑐𝜏^𝛼Γ()𝑢3−𝛼^𝑘𝑖−12^𝑛_𝑗=−𝑛𝑞_𝑗𝑖𝑢^𝑘_𝑗=2𝑏₀−𝑏₁𝑐𝜏^𝛼Γ()𝑢3−𝛼^𝑖𝑘−1−1𝑐𝜏^𝛼Γ()×3−𝛼_{𝑘−1𝑙=2}𝑏_𝑙−2−2𝑏_𝑙−1+𝑏_𝑙𝑢_𝑖^𝑘−𝑙+𝑏_𝑘−2−𝑏_𝑘−1𝑢⁰𝑥_𝑖−𝜏𝑏_𝑘−1𝜓𝑥_𝑖+12^𝑛_𝑗=−𝑛𝑞_𝑗𝑖𝑢_𝑗^𝑘−1+1𝐾𝑓^𝑘−(1/2)𝑥_𝑖, (44) where 𝑖=−𝑛,−𝑛+1,…,𝑛. To obtain a matrix representation of the above equation, we let 𝑏𝐴=⁰𝑐𝜏^𝛼Γ()𝐼3−𝛼^{(2𝑛+1)×(2𝑛+1)},1𝐵=2𝑞_{𝑗𝑖(2𝑛+1)×(2𝑛+1)},𝐶⁽¹⁾=2𝑏₀−𝑏₁𝑐𝜏^𝛼Γ()𝐼3−𝛼^{(2𝑛+1)×(2𝑛+1)},𝐶^(𝑙)𝑏=−^𝑙−2−2𝑏_𝑙−1+𝑏_𝑙𝑐𝜏^𝛼Γ()𝐼3−𝛼^{(2𝑛+1)×(2𝑛+1)},𝐶^(𝑘)𝑏=−^𝑘−2−𝑏_𝑘−1𝑐𝜏^𝛼Γ()𝐼3−𝛼^{(2𝑛+1)×(2𝑛+1)},𝑏𝐷=^𝑘−1𝑐𝜏^𝛼−1Γ()𝐼3−𝛼^{(2𝑛+1)×(2𝑛+1)},𝕌^𝑘=𝑢^𝑘_−𝑛,𝑢^𝑘_−𝑛+1,…,𝑢^𝑘_𝑛^𝑇,Ψ=𝜓𝑥_−𝑛,𝜓𝑥_−𝑛+1,…,𝜓𝑥_𝑛^𝑇,𝐹^𝑘−(1/2)=1𝐾𝑓^𝑘−(1/2)(𝑥_−𝑛),𝑓^𝑘−(1/2)(𝑥_−𝑛+1),…,𝑓^𝑘−(1/2)(𝑥_𝑛)^𝑇, (45) 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 𝑄𝕌^𝑘=𝑃, (46) where 𝑄=𝐴−𝐵,𝑃=𝐵𝕌^𝑘−1+^𝑘_𝑙=1𝐶^(𝑙)𝕌^𝑘−𝑙+𝐷Ψ+𝐹^𝑘−(1/2). (47) For the initial condition (39), we have 𝕌⁰=𝜑𝑥_−𝑛,𝜑𝑥_−𝑛+1,…,𝜑𝑥_𝑛^𝑇. (48) Consequently (46) can be solved easily for the unknown coefficients 𝕌^𝑘. Hence the approximation solution given in (41) can be obtained.

4. 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.

Lemma 6 (see [18]). For any 𝑄=𝑄₁,𝑄₂,𝑄₃,… and 𝜐, here 𝑄_𝑖,𝜐∈𝐿²(𝑎,𝑏), it can be verified that ^𝑘_𝑖=1𝑏₀𝑄_𝑖−^𝑖−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑄_𝑖−𝑙−𝑏_𝑖−1𝜐,𝑄_𝑖≥12(2−𝛼)𝑘^𝑘1−𝛼_𝑖=1𝑄_𝑖²−12𝑘^2−𝛼‖𝜐‖²,𝑘=1,2,3,…, (49) where 𝑏_𝑙 is defined in (23).

Lemma 7. Suppose that 𝑢^𝑘_𝑁, 𝑘=1,2,3,…, satisfy 𝜏^1−𝛼()𝑏𝑐Γ3−𝛼⁰𝛿_𝑡𝑢_𝑁^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑢_𝑁^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜐=𝑢_𝑁^𝑘−(1/2)_𝑥𝑥+𝑓^𝑘−(1/2),𝑢⁰_𝑁=𝜑,0<𝑥<1,𝑢^𝑘_𝑁(𝑎)=0,𝑢^𝑘_𝑁(𝑏)=0. (50)(51)(52) Then one has the estimate (𝑢^𝑘_𝑁)_𝑥²≤(𝑢⁰_𝑁)_𝑥²+𝑡_𝑘^2−𝛼()‖𝜐‖𝑐Γ3−𝛼²+𝑐Γ(2−𝛼)𝑡_𝑘^𝛼−1𝜏^𝑘_𝑙=1𝑓^{𝑘−(1/2)2}. (53)

Proof. Based on (50), we have 𝜏^1−𝛼𝑐Γ(3−𝛼)^𝑚_𝑘=1𝑏₀𝛿_𝑡𝑢_𝑁^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑢_𝑁^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜐,𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2)=^𝑚_𝑘=1𝑢_𝑁^𝑘−(1/2)_𝑥𝑥,𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2)+^𝑚_𝑘=1𝑓^𝑘−(1/2),𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2). (54) Applying Lemma 6 to the term in left side of (54), we obtain 𝜏^1−𝛼𝑐Γ(3−𝛼)^𝑚_𝑘=1𝑏₀𝛿_𝑡𝑢_𝑁^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑢_𝑁^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜐,𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2)≥1()𝑡2𝑐Γ2−𝛼^𝑚1−𝛼𝜏^𝑚_𝑘=1𝛿_𝑡𝑢_𝑁^{𝑘−(1/2)2}−1()𝑡2𝑐Γ3−𝛼^𝑚2−𝛼‖𝜐‖². (55) On the other side, based on the fact that 𝛿_𝑡𝑢_𝑁^𝑘−(1/2)(𝑎)=𝛿_𝑡𝑢_𝑁^𝑘−(1/2)(𝑏)=0, a straightforward calculation of the right terms in (54) gives ^𝑚_𝑘=1𝑢_𝑁^𝑘−(1/2)_𝑥𝑥,𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2)=−𝜏^𝑚_𝑘=1𝑢_𝑁^𝑘−(1/2)_𝑥,𝛿_𝑡𝑢_𝑁^𝑘−(1/2)_𝑥1=−2^𝑚_𝑘=1𝑢^𝑘_𝑁𝑥+𝑢_𝑁^𝑘−1_𝑥,𝑢^𝑘_𝑁𝑥−𝑢_𝑁^𝑘−1_𝑥1=−2𝑢^𝑚_𝑁𝑥²−𝑢⁰_𝑁𝑥²._𝑚𝑘=1𝑓^𝑘−(1/2),𝜏𝛿_𝑡𝑢_𝑁^𝑘−(1/2)≤12𝜏^𝑚_𝑘=1𝑐Γ(2−𝛼)𝑡_𝑚^𝛼−1𝑓^{𝑘−(1/2)2}+1()𝑡𝑐Γ2−𝛼^𝑚1−𝛼𝛿_𝑡𝑢_𝑁^{𝑘−(1/2)2}. (56) Substituting (55) and (56) into (54) yields inequality (53).

Suppose that 𝑢^𝑘, 𝑘=0,1,2,…, satisfy the conditions of Theorem 4. Substituting the approximation solution 𝑢^𝑘_𝑁 given in (41) into (34) leads to 𝜏^1−𝛼()𝑏𝑐Γ3−𝛼⁰𝛿_𝑡𝑢_𝑁^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑢_𝑁^{𝑘−𝑙−(1/2)}−𝑏_𝑘−1𝜓=𝑢_𝑁^𝑘−(1/2)_𝑥𝑥+1𝐾𝑓^𝑘−(1/2)+𝑅^𝑘_𝑁,𝑘=1,2,3,…. (57) Subtracting (57) from (34), we obtain 𝜏^1−𝛼()×𝑏𝑐Γ3−𝛼₀1𝜏𝑢^𝑘−𝑢^𝑘_𝑁−𝑢^𝑘−1−𝑢_𝑁^𝑘−1−1𝜏^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑢^𝑘−𝑙−𝑢_𝑁^𝑘−𝑙−𝑢^{𝑘−𝑙−1}−𝑢_𝑁^{𝑘−𝑙−1}=𝑢^𝑘−(1/2)_𝑥𝑥−𝑢_𝑁^𝑘−(1/2)_𝑥𝑥−𝑅^𝑘_𝑁,𝑘=1,2,3,…. (58) Based on Theorem 4, we have 𝑅^𝑘_𝑁1=−𝑐Γ(3−𝛼)𝜏^𝛼√𝑛𝑒^{−𝜋𝑑𝛽𝑛}×𝑏₀𝑑_𝑘−𝑑_𝑘−1−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝑑_𝑘−𝑙−𝑑_{𝑘−𝑙−1}+𝑑^𝑘+1𝑛^(3/2)𝑒^{−𝜋𝑑𝛽𝑛}, (59) where 𝑑₀,𝑑₁,…,𝑑_𝑘+1 are constants.

Theorem 8. Let 𝑢^𝑘_𝑁, 𝑘=0,1,2,…, be the approximation solution given in (41) and let 𝑢^𝑘 satisfy the conditions of Theorem 4. Then one has 𝑢^𝑘_𝑁²≤𝐶(𝑢⁰_𝑁)_𝑥²+𝑡_𝑘^2−𝛼()𝜓𝑐Γ3−𝛼²+𝑐Γ(2−𝛼)𝑡_𝑘^𝛼−1𝜏×^𝑘_𝑙=11𝐾𝑓^𝑙−(1/2)+𝑅_𝑁^{𝑙−(1/2)2}, (60) where 𝐶 is constant.

Proof. Considering (57) and noting that 𝑢⁰_𝑁=𝜑 and 𝑢^𝑘_𝑁(𝑎)=𝑢^𝑘_𝑁(𝑏)=0, 𝑘=0,1,2,…, we can easily obtain this result from Lemma 7 and the Poincaré inequality.

Theorem 9. Let the problem (19)–(21) have the exact solution 𝑈^𝑘 and let 𝑢^𝑘_𝑁 be the approximation solution given in (41), where 𝑘=0,1,2,…. If 𝑢^𝑘 satisfy the conditions of Theorem 4, then for 𝑘𝜏≤𝑇, one has 𝑒^𝑘≤Γ(2−𝛼)𝑇^𝛼(𝑏−𝑎)𝑐₁𝜏^3−𝛼+𝑐₂𝜏^−𝛼𝑛^(3/2)𝑒^{−𝜋𝑑𝛽𝑛}, (61) where 𝑒^𝑘=𝑈^𝑘−𝑢^𝑘_𝑁.

Proof. Subtracting (57) from (32), we obtain the error equation 𝜏^1−𝛼()𝑏𝑐Γ3−𝛼⁰𝛿_𝑡𝑒^𝑘−(1/2)−^𝑘−1_𝑙=1𝑏_𝑙−1−𝑏_𝑙𝛿_𝑡𝑒^{𝑘−𝑙−(1/2)}=𝑒^𝑘−(1/2)_𝑥𝑥+𝑅^𝑘_𝑇−𝑅^𝑘_𝑁, (62) where 𝑘=1,2,3,…. Noting that 𝑒⁰=0 and 𝑒^𝑘(𝑎)=𝑒^𝑘(𝑏)=0, 𝑘=1,2,…, and applying Lemma 7, we have 𝑒^𝑘_𝑥²≤𝑐Γ(2−𝛼)𝑡_𝑘^𝛼−1𝜏^𝑘_𝑙=1𝑅^𝑙_𝑇−𝑅^𝑙_𝑁². (63) Based on (33) and (59), we obtain 𝑒^𝑘_𝑥≤𝑐₃Γ(2−𝛼)𝑇^𝛼()×𝜏𝑏−𝑎_3−𝛼+𝑐₄𝜏^−𝛼𝑛^(3/2)𝑒^{−𝜋𝑑𝛽𝑛}. (64) Applying the Poincaré inequality yields the needed result.

Remark 10. In Theorem 9, the error estimate formula contains 𝑐₂𝜏^−𝛼𝑛^(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.

5. 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: 𝜕^𝛼𝑢(𝑥,𝑡)𝜕𝑡^𝛼=𝜕²𝑢(𝑥,𝑡)𝜕𝑥²()𝑢()()+sin𝜋𝑥,0<𝑥<1,0<𝑡≤1,𝑥,0=0,𝜕𝑢𝑥,0𝑢()()𝜕𝑡=0,0<𝑥<1,0,𝑡=0,𝑢1,𝑡=0,0<𝑡≤1. (65) The exact solution of the above problem is [33] 𝑢()=1𝑥,𝑡𝜋²1−𝐸_𝛼−𝜋²𝑡^𝛼sin(𝜋𝑥), (66) 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 𝐿²-errors and 𝐿^∞-errors. Figure 1 gives the 3D 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.

(a)

(b)

Figure 1 

The numerical solutions with 𝛼=1.3 (a) and 𝛼=1.7 (b).

(a)

(b)

Figure 2 

Comparison of the numerical and exact solutions for several fixed time instants with 𝛼=1.3 (a) and 𝛼=1.7 (b).

(a)

(b)

Figure 3 

Plot of the absolute errors with 𝛼=1.3 (a) and 𝛼=1.7 (b).

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 𝐿². All the numerical results reported in Figures 4 and 5, and Tables 1 and 2 are evaluated at 𝑡=0.8.

Table 1 

Temporal errors and convergence orders with 𝛼=1.3.

Table 2 

Temporal errors and convergence orders with 𝛼=1.7.

(a)

(b)

Figure 4 

Temporal errors in the 𝐿∞ and 𝐿2 norms versus 𝜏 for several 𝛼.

(a)

(b)

Figure 5 

Spatial errors in the 𝐿∞ and 𝐿2 norms versus 𝑛 for several 𝛼.

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 𝐿² 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.

Now we check the spatial accuracy with respect to the parameter 𝑛 about the number of collocation points. For the similar reason mentioned above, we fix the time step 𝜏 sufficiently small to avoid contamination of the temporal error. In Figure 5, we present the spatial errors as a function of 𝑛 with 𝜏=0.0001, where a logarithmic scale is now used for the spatial errors axis. It is clearly observed that the spatial errors appear in an exponential decay. In this semilog representations, we observe that the error variations are essentially linear versus 𝑛 [36].

6. Conclusion

In this paper, we have developed and analyzed an efficient numerical method for a class of the initial-boundary-value problems of fractional diffusion-wave equations. Based on a finite difference scheme in time and a global sinc-collocation method in space, the problem is reduced to the solution of a system of linear algebraic equations at each time step. The proposed method leads to 3−𝛼 order accuracy in time and exponential convergence in space. The theoretical results are perfectly confirmed by the numerical experiments.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The present research is supported by the National Natural Science Foundation of China (Grant nos. 11271311 and 11371016), the Chinese Program for Changjiang Scholars and Innovative Research Team in University (PCSIRT) (Grant no. IRT1179), and the Hunan Province Innovation Foundation for Postgraduate (Grant no. CX2013B252). The authors are grateful to the anonymous reviewers for their careful comments and valuable suggestions on this paper.