A HIGH ORDER DIFFERENCE METHOD FOR FRACTIONAL SUB-DIFFUSION EQUATIONS WITH THE SPATIALLY VARIABLE COEFFICIENTS UNDER PERIODIC BOUNDARY CONDITIONS∗

In this paper, we propose a difference scheme with global convergence order O(τ + h) for a class of the Caputo fractional equation. The difficulty caused by the spatially variable coefficients is successfully handled. The unique solvability, stability and convergence of the finite difference scheme are proved by use of the Fourier method. The obtained theoretical results are supported by numerical experiments.


Introduction
Growing interests are drawn to the study of fractional differential equations (FDEs) in recent years. FDEs have been applied to many fields, such as chemical problems and electromagnetism [19], finance [22], physics [32] and so on. There are several ways to solve the FDEs, including Green function method, Laplace and Fourier transform method etc. However, most of them are still not resolved, it is thus necessary to develop numerical methods for solving FDEs.
A fractional sub-diffusion equation is an integro-partial differential equation obtained from the classical diffusion equation by replacing the first-order time derivative by a fractional derivative of order between zero and one. We refer papers [5,11,17] to get more details. One popular approximation of the time Caputo fractional derivative is the L 1 formula, see [24]. We can refer to [3,6,7,14,21,28] for more information on the application of the L 1 approximation. Based on the Grünwald-Letnikov approximation, Yuste and Acedo [30] constructed an explicit difference scheme for fractional diffusion equations, and investigated the stability using the von Neumann method. By the finite difference method in time and the Legendre spectral method in space, a high order efficient scheme was established in [14]. While, in [35], Zhuang et al. studied the stability and convergence of an implicit numerical method for the anomalous sub-diffusion equation by the energy method. A compact finite difference scheme with O(τ + h 4 ) accuracy for this equation was then presented by Cui in [3]. By a transformation to the equation in [3,35], a compact scheme for the fractional sub-diffusion equation was shown to converge with O(τ 2−α + h 4 ) in [7]. A modified anomalous time fractional sub-diffusion equation with a nonlinear source term was studied in [2,16], where a finite difference scheme of first order temporal accuracy and fourth order spatial accuracy was proposed. By using the weighted and shifted Grünwald difference operator, Wang and Vong established schemes with O(τ 2 +h 4 ) convergence order in [29], then the method was applied to solve fractional diffusion-wave equations with Neumann boundary conditions. Lately, in [1], Alikhanov proposed a new numerical differentiation formula, called the L2-1 σ formula, to approximate the Caputo fractional derivative at a special point. Based on [1], Vong et al. [25] constructed a high order scheme to solve two dimensional time-space fractional differential equations, see [13,15] for more applications about the L2-1 σ formula.
However, the works mentioned above all involve sub-diffusion equations with constant coefficients which are usually simplified models for complicated ones. In practice, many applications are described by variable diffusion coefficients [4,18,26,34]. Therefore, it is important to study problems having variable coefficients. As expected, techniques used for equations of constant coefficient cannot be applied directly to the ones with variable coefficients. The matrix corresponding to the compact operator for problems with constant coefficient is symmetric, while those for the variable coefficient ones are not and are more difficult to handle. Meanwhile, in some cases, the problem must be modeled by high-order spatial derivatives, such as the modeling formation of grooves on a flat surface ( [20,23]) and the propagation of intense laser beams in a bulk medium with Kerr nonlinearity ( [12]). Some results related to the numerical solution for high order fractional equations can be found in [8][9][10]27,31,33]. The main purpose of this paper is to solve the following problem by the finite difference method: 1) which is affected by the initial and periodic boundary conditions denotes the Caputo fractional derivative of order γ ∈ (0, 1) on u(x, t), L is the period of u(x, t) with respect to the variable x, p(x), q(x), and r(x) are assumed to be smooth enough and satisfy p(x) ≥ 0, r(x) ≤ 0.
The paper is organized as follows. A finite difference scheme will be proposed in the next section. Unique solvability, stability and convergence of the proposed scheme are analyzed in Section 3 and 4. In Section 5, numerical experiments are carried out to support the theoretical results, the paper ends with a simple conclusion.

The proposed finite difference scheme
To develop a finite difference scheme for the problem (1.1)-(1.2), we let h = L M and τ = T N be the spatial and temporal step sizes respectively, where M and N are two given positive integers. For i = 0, 1, . . . , M and n = 0, 1, . . . , N , denote For any u ∈ V h , considering the periodic conditions, we have u n i = u n i±M . We then present a sequence {c [1], which will be used later. Let The discretization in temporal direction is based on the following lemma in [1].
The following lemma gives theoretical support for truncation error in spatial direction, the proof is just by the Taylor expansion, so we omit it here.
and define three operators A 1 , A 2 , A 3 as follows: , by the Taylor expansion, it holds that , based on Lemma 2.1 and Lemma 2.2, we propose the following scheme for the problem (2.1) One can easily check that our proposed scheme has truncation error equal to O(τ 2 + h 4 ).

Unique solvability of the difference scheme
In this section, we study the unique solvability for the proposed scheme (2.1) by the Fourier analysis method. We first give some basic preparations. Define v n (x) and z n−1+σ (x) be two periodic functions with period L, and denote restrictions on The corresponding Fourier series of v n (x) and z n−1+σ (x) can be given as Let us recommend the following discrete L 2 −norm by According to the Parseval identity, we can get the following equalities: We now turn to analysis of the unique solvability of the difference scheme (2.1), which can be written as We can then conclude the following results (3.1)

Note that
Ae βxj =e βxj (y 1 + y 2 j),   To prove the unique solvability of finite difference scheme (2.1), we only need to prove that the homogeneous difference equation has only zero solution: (3.5) Since p i ≥ 0, r i ≤ 0, we can get y 1 ≤ 0. By a similar procedure to derive (3.4), the following formula can be obtained:

Stability and convergence
In this part, we will study the stability and convergence of the finite difference scheme (2.1).

Lemma 4.1. For equation (3.4), it holds that
Proof. We can prove (4.1) by the induction. When n = 1, due to (3.4), we can conclude that Noticing that y 1 ≤ 0, σ > 1 2 and Now we turn to the general n. Suppose that (4.1) holds for 1, 2, . . . , n − 1, then Due to The desired result is proved. We now conclude the following stability and convergence result, respectively.
Proof. Based on Lemma 4.1, and the Minkowski inequality, we can deduce , the conclusion then follows by Parseval identity.
Then there exists a positive constant C such that Proof. We can easily get the following error equation . Based on the previous analysis and Theorem 4.1, we can obtain

Numerical experiments
In this section, we carry out numerical experiments on the proposed finite difference scheme to illustrate our theoretical statements. All our tests were done in MATLAB 2014b. The L 2 norm errors between the exact and the numerical solutions are shown in the following tables. Furthermore, the temporal convergence order, denoted by for sufficiently small h, the spatial convergence order, denoted by when τ is sufficiently small, are reported.
Example 5.1. We consider the following problem: The exact solution for this problem is u(x, t) = (t 2 + 1) sin(2πx), x ∈ R. The convergence order in temporal direction with h = 1 200 is reported in Table 1, while in Table 2, the convergence order in spatial direction with τ = 1 10000 , γ = 0.5 is listed. The convergence orders of the numerical results match that of the theoretical ones.

Conclusions
In this paper, the numerical solution for a fourth-order fractional sub-diffusion equations with the spatially variable coefficients under periodic boundary conditions is considered, where the global convergence order O(τ 2 + h 4 ) is obtained. The difficulty caused by the spatially variable coefficients and high order spatial derivatives is carefully handled. The unique solvability, stability and convergence of the finite difference scheme are proved by the Fourier method. Numerical experiments are carried out to justify the theoretical result.