PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 99(113) (2016), 1–13 DOI: 10.2298/PIM1613001D

FACTORIZED DIFFERENCE SCHEME FOR TWO-DIMENSIONAL SUBDIFFUSION EQUATION IN NONHOMOGENEOUS MEDIA

Aleksandra Delić, Sandra Hodžić, and Boško S. Jovanović

Abstract. A factorized finite-difference scheme for numerical approximation of initial-boundary value problem for two-dimensional subdiffusion equation in nonhomogeneous media is proposed. Its stability and convergence are investigated. The corresponding error bounds are obtained.

1. Introduction

Fractional partial differential equations have broad applications in mathematics, sciences and engineering (see [9, 15, 16]), such as anomalous transport in disordered systems, processes in viscoelastic and porous media, biological and social phenomena etc. The important characteristic of fractional differential equations is their non-local property. This means that the next state of the system depends not only upon its current state but also upon all of its previous states. Thus the fractional-order models are more realistic and it is the main reason of their popularity. On the other side, this feature makes the design of accurate and fast numerical methods difficult.

In this paper we consider the first initial-boundary value problem for twodimensional fractional in time diffusion equation with variable coefficients. This equation is commonly called subdiffusion equation. The problem is approximated by factorized finite difference scheme, which belongs to the type of alternating direction implicit (ADI) schemes (see [17]). This scheme combines the advantages of explicit and implicit scheme–efficiency and stability. Analogous result for the problem with constant coefficients is obtained in [6]. In [7] the same problem is approximated by additive scheme–another type of ADI schemes. In [3, 11] multidimensional evolution equations with fractional in space derivatives are approximated using diverse ADI schemes.

2010 Mathematics Subject Classification: 35R11, 65N12, 65N15.

Key words and phrases: fractional derivative, subdiffusion, finite differences, factorized scheme, stability, convergence rate.

Supported by Ministry of Education, Science and Technological Development of Republic of Serbia, project 174015.

Communicated by Gradimir Milovanović.

The paper is organized as follows. In Section 2 we introduce Riemann–Liouville and Caputo fractional derivatives. In Section 3 we define the first initial-boundary value problem for two-dimensional fractional in time diffusion equation and prove existence and uniqueness of its solution. In Section 4 we define factorized difference scheme approximating considered problem and prove its stability. In Section 5 we investigate the convergence of the proposed factorized difference scheme.

2. Fractional derivatives

Let u be a function defined on interval [a, b] and k − 1 ⩽ σ < k, k ∈ ℕ. The left Riemann–Liouville fractional derivative of order σ is defined as

(2.1)  D_{a+}^{σ} u(t) = \frac{1}{Γ(k − σ)} \frac{d^{k}}{dt^{k}} ∫_{a}^{t} \frac{u(s)}{(t − s)^{σ+1−k}} ds, t ⩾ a,

where the Γ(·) denotes the Gamma function. The right Riemann–Liouville fractional derivative is defined analogously

D_{b−}^{σ} u(t) = \frac{(−1)^{k}}{Γ(k − σ)} \frac{d^{k}}{dt^{k}} ∫_{t}^{b} \frac{u(s)}{(s − t)^{σ+1−k}} ds, t ⩽ b.

For σ = k −1 from (2.1) immediately follows that D_{a+}^{k−1} u(t) = u^{(k−1)} (t). Moreover, under some natural assumptions [16], we have lim_{σ→k} D_{a+}^{σ} u(t) = u^{(k)}(t).

Commuting the derivative and integral in definition (2.1) one obtains the socalled Caputo fractional derivative

CD_{a+}^{σ} u(t) = \frac{1}{Γ(k − σ)} ∫_{a}^{t} \frac{u^{(k)}(s)}{(t − s)^{σ+1−k}} ds.

The two definitions are linked by the following relationship

D_{a+}^{σ} u(t) = ^{C}D_{a+}^{σ} u(t) + ∑_{j=0}^{k−1} u^{(j)} (a)\frac{(x − a)^{j−σ}}{Γ(j − σ + 1)} .

In particular, D_{a+}^{σ} u(t) = ^{C}D_{a+}^{σ} u(t) if u(a) = u′(a) = ⋯= u^{(k−1)} (a) = 0.

Unlike classical derivatives, fractional derivatives satisfy the semigroup property under certain additional assumptions [16], for example for continuous functions:

(2.2)  D_{a+}^{σ} D_{a+}^{ϱ} u(t) = D_{a+}^{σ+ϱ} u(t) if 0 < σ, ϱ < 1, u(a) = 0.

Let 0 < σ < 1, and let u(t) and v(t) be continuously differentiable functions. Then, using the relationship between the Riemann–Liouville and Caputo fractional derivatives one easily obtains

(2.3)  (D_{a+}^{σ} u, v)_{L^{2} (a,b)} = (u, D_{b−}^{σ} v)_{L^{2} (a,b)} .

Let us mention another result that will be used in the sequel. Let σ > 0 and let u be infinitely differentiable function in ℝ, with suppu ⊂ (a, b). Then (see [5]):

(2.4)  (D_{a+}^{σ} u, D_{b−}^{σ} u)_{L^{2} (a,b)} = cosπσ∥D_{a+}^{σ} u∥_{L^{2} (a,+∞)}^{2} .

For the functions of many variables, the partial fractional derivatives are defined in an analogous manner, for example

D_{t,a+}^{σ} u(x, t) = \frac{1}{Γ(k − σ)} \frac{∂^{k}}{∂t^{k}} ∫_{a}^{t} \frac{u(x, s)}{(t − s)^{σ+1−k}} ds, k − 1 < σ < k, k ∈ ℕ.

3. Problem formulation

We shall consider the time fractional diffusion equation

(3.1)  D_{t,0+}^{α} u + Lu = f(x, t), 0 < α < 1, x = (x_{1}, x_{2}) ∈ Ω, t ∈ (0, T),

where Ω = (0, 1) × (0, 1), Q = Ω × (0, T) and

Lu = − ∑_{i,j=1}^{2} \frac{∂}{∂x_{i}} ( a_{ij} \frac{∂u}{∂x_{j}} ) + au,

are subject to homogeneous boundary and initial conditions

(3.2)  u(x, t) = 0,  x ∈ ∂Ω,  t ∈ (0, T),
(3.3)  u(x, 0) = 0,  x ∈ Ω̄.

We assume that the coefficients of the differential operator Lu satisfy the standard ellipticity assumptions

(3.4)  \begin{matrix}a_{ij} , a ∈ L^{∞} (Ω), & a ⩾ 0, & a_{ij} = a_{ji} , \\ ∑_{i,j=1}^{2} a_{ij} ξ_{i}ξ_{j} ⩾ c_{0} ∑_{i=1}^{2} ξ_{i}^{2}, & x ∈ Ω, & ξ = (ξ_{1}, ξ_{2}) ∈ ℝ^{2}, c_{0} > 0. \end{matrix}

Initial-boundary value problem (3.1)–(3.3) is often called sub-diffusion problem.

With C^{k}(Ω) and C^{k}(Ω̄) we denote the spaces of k-fold differentiable functions in Ω and Ω̄, respectively. In particular, Ċ^{∞} (Ω) = C_{0}^{∞} (Ω) stand for the space of infinitely differentiable functions with compact support in Ω. As usual, the space of measurable functions whose squares are Lebesgue integrable in Ω is denoted by L^{2}(Ω). We also use W^{σ,p} (Ω) to denote the Sobolev spaces [1]. By Ẇ^{σ,p}(Ω) = W_{0}^{σ,p} (Ω) we denote the closure of Ċ^{∞} (Ω) with the respect to the norm of W^{σ,p} (Ω). In particular, for p = 2 we set H^{σ}(Ω) = W^{σ,2} (Ω) and Ḣ^{σ} (Ω) = Ẇ^{σ,2}(Ω).

For σ > 0 we set

|u|_{C_{+}^{σ} [a,b]} = ∥D_{a+}^{σ} u∥_{C[a,b]}, |u|_{C_{−}^{σ} [a,b]} = ∥D_{b−}^{σ} u∥_{C[a,b]},
∥u∥_{C_{±}^{σ} [a,b]} = (∥u∥_{C^{[σ]^{−}} [a,b]}^{2} + |u|_{C_{±}^{σ} [a,b]}^{2})^{1/2} ,
|u|_{H_{+}^{σ} (a,b)} = ∥D_{a+}^{σ} u∥_{L^{2} (a,b)} , |u|_{H_{−}^{σ} (a,b)} = ∥D_{b−}^{σ} u∥_{L^{2} (a,b)} ,
∥u∥_{H_{±}^{σ} (a,b)} = (∥u∥_{H^{[σ]^{−}} (a,b)}^{2} + |u|_{H_{±}^{σ} (a,b)}^{2} )^{1/2} ,

where [σ]^{−} denotes the largest integer < σ. Then we define C_{±}^{σ} [a, b] as the space of functions u ∈ C^{[σ]^{−}} [a, b] with the finite norms ∥u∥_{C_{±}^{σ} [a,b]}. The space H_{±}^{σ} (a, b) is defined analogously, while the space Ḣ_{±}^{σ} (a, b) is defined as the closure of Ċ^{∞} (a, b) with respect to the norm ∥ · ∥_{H_{±}^{σ} (a,b)} . Because for σ = k ∈ ℕ ∪ {0} fractional derivative reduces to standard k-th derivative, we have C_{±}^{k} [a, b] = C^{k}[a, b] and H_{±}^{k} (a, b) = H^{k}(a, b).

The following result holds:

Lemma 3.1. (See [12]) For σ > 0, σ ≠ k + 1/2, k ∈ ℕ, the spaces Ḣ_{+}^{σ} (a, b), Ḣ_{−}^{σ} (a, b) and Ḣ^{σ} (a, b) are equal and their norms are equivalent.

For the vector valued functions mapping real interval [0, T] or (0, T) into Banach space X we introduce the spaces C^{k}([0, T], X), k ∈ ℕ ∪ {0} and H^{σ}((0, T), X), σ ⩾ 0, in the usual way [13]. In analogous manner we define the spaces C_{±}^{σ} ([0, T], X) and H_{±}^{σ} ((0, T), X).

Taking the inner product of equation (3.1) with test function v and using partial integration and relations (2.2) and (2.3) one obtains the following weak formulation of the problem (3.1)–(3.3): find u ∈ Ḣ^{1,α/2} (Q) = L^{2}((0, T), Ḣ^{1}(Ω)) ∩ Ḣ^{α/2} ((0, T), L^{2}(Ω)) such that

a(u, v) = l(v), ∀v ∈ Ḣ^{1,α/2} (Q),

where

a(u, v) = (D_{t,0+}^{α/2} u, D_{t,T −}^{α/2} v)_{L^{2} (Q)} + ∑_{i,j=1}^{2} ( a_{ij} \frac{∂u}{∂x_{j}} , \frac{∂v}{∂x_{i}} )_{L^{2} (Q)} + (au, v)_{L^{2} (Q)}

and

l(v) = (f, v)_{L^{2} (Q)} .

Theorem 3.1. Let α ∈ (0, 1), f ∈ L^{2}(Q) and let the assumptions (3.4) hold. Then the problem (3.1)–(3.3) is well posed in Ḣ^{1,α/2} (Q) and its weak solution satisfies a priori estimate

∥u∥_{H^{1,α/2} (Q)} ⩽ C∥f∥_{L^{2} (Q)} .

The proof follows immediately using relation (2.4) and the Lax–Milgram lemma.

Remark 3.1. From Theorem 3.1 immediately follows a priori estimate

(3.5)  ∥u∥_{B^{1,α/2} (Q)} ⩽ C∥f∥_{L^{2} (Q)}

in weaker norm [13]

∥u∥_{B^{1,α/2} (Q)}^{2} = ∫_{0}^{T} [(T − t)^{−α} ∥u(·, t)∥_{L^{2} (Ω)}^{2} + ∥u(·, t)∥_{H^{1} (Ω)}^{2} ]dt.

Theorem 3.2. Let the assumptions of Theorem 3.1 hold and let a_{ij} ∈ W^{1,∞} (Ω). Then solution of the problem (3.1)–(3.3) belongs to the space H_{+}^{2,α}(Q) ∩ Ḣ^{1,α/2} (Q) and the a priori estimate holds

∥u∥_{H_{+}^{2,α} (Q)} ⩽ C∥f∥_{L^{2} (Q)} .

The proof follows taking the inner product of (3.1) with Lu, using relation (2.4) and the so-called ‘second fundamental inequality’ [10]

∥u∥_{H^{2} (Ω)}^{2} ⩽ C(∥Lu∥_{L^{2} (Ω)}^{2} + ∥u∥_{H^{1} (Ω)}^{2} ).

In the sequel, we shall assume that a_{ij} , a ∈ C(Ω̄) instead of L^{∞} (Ω). Through the paper by C we shall denote the positive generic constant which may take different values in different formulas.

4. Finite difference approximation

In the area Q̄ = [0, 1]×[0, 1]×[0, T], we define the uniform mesh Q̄_{hτ} = ω̄_{h} ×ω̄_{τ} , where ω̄_{h} = {x = (n_{1}h, n_{2}h) | n_{1}, n_{2} = 0, 1, ..., N; h = 1/N} and ω̄_{τ} = {t = t_{m} = mτ | m = 0, 1, …, M; τ = T/M}. We also define ω_{h} = ω̄_{h} ∩ Ω, γ_{h} = ω̄_{h} ∖ ω_{h}, ω_{0h} = ω̄_{h} ∩ (0, 1] × (0, 1], ω_{1h} = ω̄_{h} ∩ (0, 1] × (0, 1), ω_{2h} = ω̄_{h} ∩ (0, 1) × (0, 1], ω_{τ} = ω̄_{τ} ∩ (0, T), ω_{τ}^{−} = ω̄_{τ} ∩ [0, T) and ω_{τ}^{+} = ω̄_{τ} ∩ (0, T]. We shall use the standard notation from the theory of the finite difference schemes (see [17]):

v = v(x, t), v̂ = v(x, t + τ), v̌ = v(x, t − τ), v^{m} = v(x, t_{m}), x ∈ ω̄_{h} ,

v_{x_{i}} = \frac{v(x + he_{i}, t) − v(x, t)}{h} = vx̄_{i} (x + he_{i}, t), i = 1, 2,

v_{t} = \frac{v(x, t + τ) − v(x, t)}{τ} = vt̄(x, t + τ) = v̂t̄,

where e_{i} denotes the unit vector of the axis 0x_{i}.

For a function u defined on Q̄ which satisfies zero initial condition, we approximate the left Riemann–Liouville fractional derivative D_{t,0+}^{α} u(x, t_{m}), 0 < α < 1, by (see [4]):

(D_{t,0+,τ}^{α} u)^{m} = \frac{1}{Γ(2 − α)} ∑_{l=0}^{m−1} (t_{m−l}^{1−α} − t_{m−l−1}^{1−α} )u_{t}^{l}.

The following result holds:

Lemma 4.1. (See [18]) Suppose that α ∈ (0, 1), u ∈ C^{2}([0, t], C(Ω̄)) and t ∈ ω_{τ}^{+} . Then

|D_{t,0+}^{α} u − D_{t,0+,τ}^{α} u| ⩽ \frac{τ^{2−α}}{1 − α} [ \frac{1 − α}{12} + \frac{2^{2−α}}{2 − α} − (1 + 2^{−α} )] max_{Q̄_{t}} | \frac{∂^{2}u}{∂t^{2}} |,

where denoted Q_{t} = (0, t) × Ω.

We approximate initial-boundary value problem (3.1)–(3.3) with the following factorized finite difference scheme:

(4.1)  ((I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} v)^{m} + L_{h}v^{m−1} = f̄^{m} , x ∈ ω_{h},

m = 1, 2, …, M, subject to zero boundary and initial conditions:

(4.2)  v(x, t) = 0, (x, t) ∈ γ_{h} × ω_{τ}^{+} ,
(4.3)  v(x, 0) = 0, x ∈ ω̄_{h} ,

where I is the identity operator, θ is positive parameter, A_{i}v = −v_{x_{i} x̄_{i}} , i = 1, 2, and

L_{h}v = −\frac{1}{2} ∑_{i,j=1}^{2} [(a_{ij} v_{x_{j}} )x̄_{i} + (a_{ij} vx̄_{j} )_{x_{i}} ] + av.

When the right-hand side f is a continuous function, we set f̄ = f, otherwise we must use some averaged value, for example f̄ = T_{1}T_{2}f, where T_{i} are Steklov averaging operators:

T_{i}f(x, t) = ∫_{−1/2}^{1/2} f(x + hse_{i}, t) ds, i = 1, 2.

Let us note, that the finite difference scheme (4.1)–(4.3) is numerically efficient, unlike the standard implicit scheme [19]. Indeed, to compute the values of solution v at the time layer t = t_{m} it is necessary to invert the operators (I + θτ^{α} A_{1}) and (I + θτ^{α} A_{2}) . By suitable ordering of mesh nodes in ω_{h} each of these operators can be represented by a tridiagonal matrix. In this way, the required values of solution are obtained by two applications of the Thomas algorithm.

We define the following discrete inner products and norms:

(v, w)_{h} = (v, w)_{L^{2} (ω_{h} )} = h^{2} ∑_{x∈ω_{h}} vw, ∥v∥_{h} = ∥v∥_{L^{2} (ω_{h} )} = (v, v)_{h}^{1/2} ,

(v, w)_{ih} = (v, w)_{L^{2} (ω_{ih} )} = h^{2} ∑_{x∈ω_{ih}} vw, ∥v∥_{ih} = ∥v∥_{L^{2} (ω_{ih} )} = (v, v)_{ih}^{1/2} , i = 0, 1, 2,

|v|_{H^{1} (ω_{h} )}^{2} = ∑_{i=1}^{2} ∥vx̄_{i} ∥_{ih}^{2} , ∥v∥_{H^{1} (ω_{h} )}^{2} = |v|_{H^{1} (ω_{h} )}^{2} + ∥v∥_{h}^{2},

|v|_{H^{2} (ω_{h} )}^{2} = ∑_{i=1}^{2} ∥v_{x̄_{i} x_{i}} ∥_{h}^{2} + 2∥v_{x̄_{1} x̄_{2}} ∥_{0h}^{2}, ∥v∥_{H^{2} (ω_{h} )}^{2} = |v|_{H^{2} (ω_{h} )}^{2} + ∥v∥_{H^{1} (ω_{h} )}^{2} ,

∥v∥_{L^{2} (Q_{hτ} )}^{2} = τ ∑_{m=1}^{M} ∥v^{m}∥_{h}^{2}, ∥v∥_{L^{2} (Q_{ihτ} )}^{2} = τ ∑_{m=1}^{M} ∥v^{m}∥_{ih}^{2} , i = 0, 1, 2,

∥v∥_{B^{1,α/2} (Q_{hτ} )}^{2} = τ ∑_{m=1}^{M} [(D_{t,0+,τ}^{α} (∥v∥_{h}^{2}))^{m} + ∥v^{m}∥_{H^{1} (ω_{h} )}^{2} ],

∥v∥_{H_{+}^{2,α} (Q_{hτ} )}^{2} = τ ∑_{m=1}^{M} [∥(D_{t,0+,τ}^{α} v)^{m}∥_{h}^{2} + ∥v^{m}∥_{H^{2} (ω_{h} )}^{2} ].

For α ∈ (0, 1) and every function v(t) defined on the mesh ω̄_{τ} , which satisfies the initial condition v(0) = 0, the following equality is valid (see [4])

(4.4)  τ ∑_{m=1}^{M} (D_{t,0+,τ}^{α} (v^{2}))^{m} = \frac{1}{Γ(2 − α)} ∑_{m=1}^{M} (t_{M−m+1}^{1−α} − t_{M−m}^{1−α} ) (v^{m})^{2}.

In particular, from here follows that the norm ∥v∥_{B^{1,α/2} (Q_{hτ} )} is well defined.

Lemma 4.2. (See [7]) For 0 < α < 1 and any function v(t) defined on the mesh ω̄_{τ} the following inequality is valid

(4.5)  v^{m}(D_{t,0+,τ}^{α} v)^{m} ⩾ \frac{1}{2}(D_{t,0+,τ}^{α} (v^{2}))^{m} + \frac{τ^{2−α} (1 − 2^{−α} )}{Γ(2 − α)} (v_{t}^{m−1} )^{2},

for m = 1, 2, …, M.

Theorem 4.1. Let α ∈ (0, 1) and θ ⩾ \frac{Γ(2−α)}{1−2^{−α}} max_{ij} ∥a_{ij} ∥_{C(Ω̄ )}. Then for sufficiently small τ finite difference scheme (4.1)–(4.3) is absolutely stable and its solution satisfies the following a priori estimate:

(4.6)  ∥v∥_{B^{1,α/2} (Q_{hτ} )} ⩽ C∥f̄∥_{L^{2} (Q_{hτ} )} .

Proof. Taking the inner product of (4.1) with v^{m}, we obtain

(v^{m}, BD_{t,0+,τ}^{α} v^{m})_{h} + (v^{m}, L_{h}v^{m})_{h} − (v^{m}, L_{h}(v^{m} − v^{m−1} ))_{h} = (v^{m}, f̄^{m} )_{h},

where denoted B = (I + θτ^{α} A_{1})(I + θτ^{α} A_{2}). The operators L_{h} and B are positive and selfadjoint, so the corresponding energy norms (see [17]) ∥v∥_{L_{h}} = (L_{h}v, v)_{h}^{1/2} and ∥v∥_{B} = (Bv, v)_{h}^{1/2} are well defined.

Using inequality (4.5), Cauchy–Schwarz and ε inequalities and taking into account that v^{m} = \frac{1}{2}(v^{m} + v^{m−1} ) + \frac{1}{2}(v^{m} − v^{m−1} ), we obtain

\frac{1}{2}D_{t,0+,τ}^{α} ∥v^{m}∥_{B}^{2} + \frac{τ^{2−α} (1 − 2^{−α} )}{Γ(2 − α)} ∥v_{t}^{m−1} ∥_{B}^{2} + \frac{1}{2}∥v^{m}∥_{L_{h}}^{2} + \frac{1}{2}∥v^{m−1} ∥_{L_{h}}^{2}

− \frac{τ^{2}}{2}∥v_{t}^{m−1} ∥_{L_{h}}^{2} ⩽ \frac{1}{4ε}∥f̄^{m} ∥_{h}^{2} + ε∥v^{m}∥_{h}^{2}, ∀ε > 0.

Further we have

∥v∥_{B}^{2} = ∥v∥_{h}^{2}+θτ^{α} ∥v∥_{A_{1} +A_{2}}^{2} +θ^{2}τ^{2α} ∥v∥_{A_{1} A_{2}}^{2} = ∥v∥_{h}^{2}+θτ^{α} |v|_{H^{1} (ω_{h} )}^{2} +θ^{2}τ^{2α} ∥v_{x̄_{1} x̄_{2}} ∥_{0h}^{2}.

Using (3.4) one obtains

c_{0}|v|_{H^{1} (ω_{h} )}^{2} ⩽ ∥v∥_{L_{h}}^{2} ⩽ c_{1}|v|_{H^{1} (ω_{h} )}^{2} + c_{2}∥v∥_{h}^{2},

where c_{1} = 2 max_{ij} ∥a_{ij} ∥_{C(Ω̄ )} and c_{2} = ∥a∥_{C(Ω̄ )}.

From the last three inequalities it follows

\frac{1}{2}D_{t,0+,τ}^{α} ∥v^{m}∥_{h}^{2} + ( θ \frac{(1 − 2^{−α} )}{Γ(2 − α)} − \frac{c_{1}}{2} ) τ^{2}|v_{t}^{m−1} |_{H^{1} (ω_{h} )}^{2}

+ ( \frac{1 − 2^{−α}}{Γ(2 − α)} − \frac{c_{2}}{2} τ^{α} ) τ^{2−α} ∥v_{t}^{m−1} ∥_{h}^{2} + \frac{c_{0}}{2} |v^{m}|_{H^{1} (ω_{h} )}^{2} + \frac{c_{0}}{2} |v^{m−1} |_{H^{1} (ω_{h} )}^{2}

⩽ \frac{1}{4ε}∥f̄^{m} ∥_{h}^{2} + ε∥v^{m}∥_{h}^{2},

whereby for

θ ⩾ θ_{0} = \frac{c_{1}Γ(2 − α)}{2(1 − 2^{−α} )} = \frac{Γ(2 − α)}{1 − 2^{−α}} max_{ij}∥a_{ij} ∥_{C(Ω̄ )}

and

τ ⩽ τ_{0} = ( \frac{2(1 − 2^{−α} )}{c_{2}Γ(2 − α)} )^{1/α} = ( \frac{2(1 − 2^{−α} )}{Γ(2 − α)∥a∥_{C(Ω̄ )}} )^{1/α}

one obtains

\frac{1}{2}D_{t,0+,τ}^{α} ∥v^{m}∥_{h}^{2} + \frac{c_{0}}{2} |v^{m}|_{H^{1} (ω_{h} )}^{2} + \frac{c_{0}}{2} |v^{m−1} |_{H^{1} (ω_{h} )}^{2} ⩽ \frac{1}{4ε}∥f̄^{m} ∥_{h}^{2} + ε∥v^{m}∥_{h}^{2}.

Finally, using the discrete Poincaré inequality (see [17])

∥v∥_{h} ⩽ \frac{1}{4} |v|_{H^{1} (ω_{h} )} ,

setting ε = 4c_{0} and summing for m = 1, 2, …, M, we obtain a priori estimate (4.6) with

C = \sqrt{\frac{1}{8c_{0} min\{1, \frac{8c_{0}}{17} \}}}  □

Note that inequality (4.6) is discrete analogue of a priori estimate (3.5). The discrete version of Theorem 3.2 we shall prove in the special case when L_{h} is finite difference approximation of Laplace operator.

Theorem 4.2. Let α ∈ (0, 1) and θ ⩾ \frac{Γ(2−α)}{2(1−2^{−α} )} . Then finite difference scheme (4.1)–(4.3) with L_{h}v = Av = (A_{1} + A_{2})v = −v_{x_{1} x̄_{1}} − v_{x_{2} x̄_{2}} is absolutely stable in the norm H_{+}^{2,α} (Q_{hτ}) and its solution satisfies a priori estimate

(4.7)  ∥v∥_{H_{+}^{2,α} (Q_{hτ} )} ⩽ C∥f̄∥_{L^{2} (Q_{hτ} )} .

Proof. Taking the inner product of (4.1) with Av^{m}, we obtain

(Av^{m}, BD_{t,0+,τ}^{α} v^{m})_{h} + (Av^{m}, Av^{m})_{h} − (Av^{m}, A(v^{m} − v^{m−1} ))_{h} = (Av^{m}, f̄^{m} )_{h},

whereby, analogously as in the proof of Lemma 4.1

\frac{1}{2}D_{t,0+,τ}^{α} ∥v^{m}∥_{AB}^{2} + \frac{τ^{2−α} (1 − 2^{−α} )}{Γ(2 − α)} ∥v_{t}^{m−1} ∥_{AB}^{2} + \frac{1}{2}∥Av^{m}∥_{h}^{2} + \frac{1}{2}∥Av^{m−1} ∥_{h}^{2}

− \frac{τ^{2}}{2}∥Av_{t}^{m−1} ∥_{h}^{2} ⩽ \frac{1}{4ε}∥f̄^{m} ∥_{h}^{2} + ε∥Av^{m}∥_{h}^{2}.

Further we have

∥v∥_{AB}^{2} ⩾ ∥v∥_{A}^{2} = |v|_{H^{1} (ω_{h} )}^{2} ,
∥v∥_{AB}^{2} ⩾ θτ^{α} ∥v∥_{A^{2}}^{2} ,
∥v∥_{A^{2}}^{2} = ∥Av∥_{h}^{2} = |v|_{H^{2} (ω_{h} )}^{2} .

From the previous inequalities it follows

\frac{1}{2}D_{t,0+,τ}^{α} |v^{m}|_{H^{1} (ω_{h} )}^{2} + ( θ \frac{(1 − 2^{−α} )}{Γ(2 − α)} − \frac{1}{2}) τ^{2}|v_{t}^{m−1} |_{H^{2} (ω_{h} )}^{2}

+ \frac{1}{2}|v^{m}|_{H^{2} (ω_{h} )}^{2} + \frac{1}{2}|v^{m−1} |_{H^{2} (ω_{h} )}^{2} ⩽ \frac{1}{4ε}∥f̄^{m} ∥_{h}^{2} + ε|v^{m}|_{H^{2} (ω_{h} )}^{2} ,

whereby for θ ⩾ \frac{Γ(2−α)}{2(1−2^{−α} )} , setting ε = 1/4 and summing for m = 1, 2, …, M, we obtain

2τ ∑_{m=1}^{M} D_{t,0+,τ}^{α} |v^{m}|_{H^{1} (ω_{h} )}^{2} + τ ∑_{m=1}^{M} |v^{m}|_{H^{2} (ω_{h} )}^{2} ⩽ 4∥f̄^{m} ∥_{L^{2} (Q_{hτ} )}^{2} .

The first sum on the left hand side is positive according to (4.4). From equation (4.1) it follows

∥(D_{t,0+,τ}^{α} v)^{m}∥_{h} ⩽ ∥Lv^{m−1} ∥_{h} + ∥f̄^{m} ∥_{h} = ∥Av^{m−1} ∥_{h} + ∥f̄^{m} ∥_{h}
= |v^{m−1} |_{H^{2} (ω_{h} )} + ∥f̄^{m} ∥_{h}.

Hence,

τ ∑_{m=1}^{M} ∥(D_{t,0+,τ}^{α} v)^{m}∥_{h}^{2} + τ ∑_{m=1}^{M} |v^{m}|_{H^{2} (ω_{h} )}^{2} ⩽ 14∥f̄^{m} ∥_{L^{2} (Q_{hτ} )}^{2}

and results follow using equivalence of seminorm |v|_{H^{2} (ω_{h} )} and norm ∥v∥_{H^{2} (ω_{h} )} for mesh functions satisfying homogeneous boundary condition v = 0 on γ_{h} (see [17]):

|v|_{H^{2} (ω_{h} )} ⩽ ∥v∥_{H^{2} (ω_{h} )} ⩽ \sqrt{1 + \frac{1}{4^{2}} + \frac{1}{16^{2}}} |v|_{H^{2} (ω_{h} )} .  □

5. Convergence of the difference scheme

Let u be solution of the initial-boundary-value problem (3.1)–(3.3) and v the solution of the difference problem (4.1)–(4.3) with f̄ = T_{1}T_{2}f. The error z = u − v is defined on the mesh ω̄_{h} × ω̄_{τ} . Putting v = −z + u into (4.1)–(4.3) we conclude that the error z satisfies the following finite difference scheme:

(5.1)  ((I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} z)^{m} + L_{h}z^{m−1} = ψ^{m},
x ∈ ω_{h}, m = 1, 2, …, M,
(5.2)  z = 0, x ∈ γ_{h}, t ∈ ω̄_{τ} ,
(5.3)  z^{0} = z(x, 0) = 0, x ∈ ω_{h},

where

ψ^{m} = (I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} u^{m} + L_{h}u^{m−1} − T_{1}T_{2}f^{m}

= ξ^{m} + η^{m} + ∑_{i,j=1}^{2} (η_{ij,x_{i}}^{m} + ζ_{ij,x_{i}}^{m} ) + ∑_{i=1}^{2} (χ_{i,x_{i}}^{m} + μ_{i,x_{i}}^{m} ),

and

ξ = D_{t,0+,τ}^{α} u − T_{1}T_{2}(D_{t,0+}^{α} u),
η = aǔ − T_{1}T_{2}(au),

η_{ij} = T_{3−i} ( a_{ij} \frac{∂u}{∂x_{j}} ) |_{(x−0.5he_{i} ,t)} − \frac{1}{2}[(a_{ij} ux̄_{j} )|_{(x,t)} + (a_{ij} u_{x_{j}} )|_{(x−he_{i} ,t)} ],

ζ_{ij} = \frac{τ}{2}[(a_{ij} u_{x̄_{j} t̄})|_{(x,t)} + (a_{ij} u_{x_{j} t̄})|_{(x−he_{i} ,t)} ],

χ_{i} = −θτ^{α} D_{t,0+,τ}^{α} ux̄_{i} ,

μ_{i} = \frac{1}{2} θ^{2}τ^{2α} D_{t,0+,τ}^{α} u_{x̄_{i} x_{3−i} x̄_{3−i}} .

Theorem 5.1. Under the assumptions of Theorem 4.1 finite difference scheme (5.1)–(5.3) is absolutely stable and the following a priori estimate holds:

(5.4) ∥z∥_{B^{1,α/2} (Q_{hτ} )} ⩽ C[ ∑_{i,j=1}^{2} (∥η_{ij} ∥_{L^{2} (Q_{ihτ} )} + ∥ζ_{ij} ∥_{L^{2} (Q_{ihτ} )} )

+ ∑_{i=1}^{2} (∥χ_{i}∥_{L^{2} (Q_{ihτ} )} + ∥μ_{i}∥_{L^{2} (Q_{ihτ} )} ) + ∥ξ∥_{L^{2} (Q_{hτ} )} + ∥η∥_{L^{2} (Q_{hτ} )} ].

The proof is analogous to the proof of Theorem 4.1, where in the estimation of truncation error terms containing finite differences the partial summation is used.

In such a way, to obtain the error bound for finite difference scheme (4.1)–(4.3) it is sufficient to estimate the right-hand side terms in (5.4).

Theorem 5.2. Let the assumptions of Theorem 4.1 hold, a_{ij} , a ∈ H^{2}(Ω) and let the solution u of initial-boundary value problem (3.1)–(3.3) belongs to the space C^{2}([0, T], C(Ω̄)) ∩ C_{+}^{α} ([0, T], H^{3}(Ω)) ∩ H^{1}((0, T), H^{2}(Ω)). Then the solution v of finite difference scheme (4.1)–(4.3) with f̄ = T_{1}T_{2}f converges to u and the following convergence rate estimate holds:

∥u − v∥_{B^{1,α/2} (Q_{hτ} )} = O(h^{2} + τ^{α}).
Proof. Let us set ξ = ξ_{1} + ξ_{2}, where
ξ_{1} = D_{t,0+,τ}^{α} u − D_{t,0+}^{α} u and ξ_{2} = D_{t,0+}^{α} u − T_{1}T_{2}D_{t,0+}^{α} u.

From Lemma 4.1 immediately it follows

∥ξ_{1}∥_{L^{2} (Q_{hτ} )} ⩽ Cτ^{2−α} ∥u∥_{C^{2} ([0,T],C(Ω̄ ))} .

From integral representation

u(x, t)−T_{1}T_{2}u(x, t) = \frac{1}{h^{2}} ∫_{x_{1} −h/2}^{x_{1} +h/2} ∫_{x_{2} −h/2}^{x_{2} +h/2} ( ∫_{x′_{1}}^{x_{1}} ∫_{x′_{2}}^{x_{2}} \frac{∂^{2}u}{∂x_{1}∂x_{2}} (x′′_{1}, x′′_{2}, t)dx′′_{2}dx′′_{1}

− ∫_{x′_{1}}^{x_{1}} ∫_{x′′_{1}}^{x_{1}} \frac{∂^{2}u}{∂x_{1}^{2}}(x′′′_{1}, x′_{2}, t)dx′′′_{1}dx′′_{1}− ∫_{x′_{2}}^{x_{2}} ∫_{x′′_{2}}^{x_{2}} \frac{∂^{2}u}{∂x_{2}^{2}}(x′_{1}, x′′′_{2}, t) dx′′′_{2}dx′′_{2}) dx′_{2}dx′_{1}

one obtains

∥ξ_{2}∥_{L^{2} (Q_{hτ} )} = ( τ ∑_{m=1}^{M} ∥ξ_{2}^{m}∥_{h}^{2}) ^{1/2} ⩽ Ch^{2}∥u∥_{C_{+}^{α} ([0,T],H^{2} (Ω))} .

Hence,

(5.5)  ∥ξ∥_{L^{2} (Q_{hτ} )} ⩽ C(τ^{2−α} ∥u∥_{C^{2} ([0,T],C(Ω̄ ))} + h^{2}∥u∥_{C_{+}^{α} ([0,T],H^{2} (Ω))} ).
Let us set η = η_{1} + η_{2}, where
η_{1} = aǔ − T_{1}T_{2}(aǔ) and η_{2} = T_{1}T_{2}(aǔ) − T_{1}T_{2}(au) = −τT_{1}T_{2}(aut̄).

The value of η_{1} at the mesh node x ∈ ω_{h} is bounded linear functional of aǔ ∈ H^{2}(e), where e = e(x) = (x_{1} − h/2, x_{1} + h/2) × (x_{2} − h/2, x_{2} + h/2), which vanishes when aǔ = 1, x_{1}, x_{2}. Using the Bramble-Hilbert lemma [2] one obtains

|η_{1}| ⩽ Ch|aǔ|_{H^{2} (e)} .

Summing this inequality over the mesh ω_{h} and using the properties of multipliers in the Sobolev spaces [14] we get

∥η_{1}∥_{h} ⩽ Ch^{2}∥aǔ∥_{H^{2} (Ω)} ⩽ Ch^{2}∥a∥_{H^{2} (Ω)} ∥ǔ∥_{H^{2} (Ω)}

and finally, after summation over the mesh ω_{τ}^{+} :

∥η_{1}∥_{L^{2} (Q_{hτ} )} ⩽ Ch^{2}∥a∥_{H^{2} (Ω)} ∥u∥_{C([0,T],H^{2} (Ω))} .

The term η_{2} can be estimated directly:

∥η_{2}∥_{L^{2} (Q_{hτ} )} ⩽ Cτ∥a∥_{C(Ω̄ )} ∥ \frac{∂u}{∂t} ∥_{L^{2} (Q)} ⩽ Cτ∥a∥_{C(Ω̄ )}∥u∥_{C^{1} ([0,T],L^{2} (Ω))} .

From the obtained inequalities it follows

(5.6) ∥η∥_{L^{2} (Q_{hτ} )} ⩽ C(h^{2}∥a∥_{H^{2} (Ω)} ∥u∥_{C([0,T],H^{2} (Ω))} + τ∥a∥_{C(Ω̄ )}∥u∥_{C^{1} ([0,T],L^{2} (Ω))} ).

Let us decompose η_{ij} in the following manner

η_{ij} = η_{ij1} + η_{ij2} + η_{ij3} + η_{ij4} , where
η_{ij1} = T_{3−i} ( a_{ij} \frac{∂u}{∂x_{j}} ) |_{(x−0.5he_{i} ,t)} − T_{3−i} (a_{ij} )|_{(x−0.5he_{i} )} T_{3−i} ( \frac{∂u}{∂x_{j}} ) |_{(x−0.5he_{i} ,t)} ,
η_{ij2} = T_{3−i} (a_{ij} ) |_{(x−0.5he_{i} )} [T_{3−i} ( \frac{∂u}{∂x_{j}} ) |_{(x−0.5he_{i} ,t)} − \frac{1}{2}(ux̄_{j} (x, t) + u_{x_{j}} (x−he_{i}, t))],
η_{ij3} = \frac{1}{2}[T_{3−i} (a_{ij} ) |_{(x−0.5he_{i} )} − \frac{1}{2}( a_{ij} (x) + a_{ij} (x − he_{i})) ]
× [ux̄_{j} (x, t) + u_{x_{j}} (x − he_{i}, t)],
η_{ij4} = −\frac{1}{4}[a_{ij} (x) − a_{ij} (x − he_{i})][ux̄_{j} (x, t) − u_{x_{j}} (x − he_{i}, t)].

The value of η_{ij1} at the mesh node x ∈ ω_{ih} is bounded bilinear functional of (a_{ij} , u) ∈ W^{1,q}(e_{i}) × W^{2,2q/(q−2)} (e_{i}), where e_{i} = e(x − 0.5he_{i}) and q > 2, which vanishes when a_{ij} = 1 or u = 1, x_{1}, x_{2}. Using the bilinear version of the Bramble– Hilbert lemma and methodology presented in [8] one obtains

|η_{ij1} (x, t)| ⩽ Ch|a_{ij} |_{W^{1,q} (e_{i} )} |u(·, t)|_{W^{2,2q/(q−2)} (e_{i} )} , q > 2.

Summing this inequality over the mesh ω_{ih} and using the Sobolev imbedding theorems [1] we get

∥η_{ij1} (·, t)∥_{ih} ⩽ Ch^{2}∥a_{ij} ∥_{W^{1,q} (Ω)} ∥u(·, t)∥_{W^{2,2q/(q−2)} (Ω)}
⩽ Ch^{2}∥a_{ij} ∥_{H^{2} (Ω)} ∥u(·, t)∥_{H^{3} (Ω)} .

Analogous inequalities hold for η_{ij2} , η_{ij3} and η_{ij4} . In such a way, after summation over the mesh ω_{τ}^{+} we obtain:

(5.7)  ∥η_{ij} ∥_{L^{2} (Q_{ihτ} )} ⩽ Ch^{2}∥a_{ij} ∥_{H^{2} (Ω)} ∥u∥_{C([0,T],H^{3} (Ω))} .

The terms ζ_{ij} , χ_{i} and μ_{i} can be estimated directly:

(5.8)  ∥ζ_{ij} ∥_{L^{2} (Q_{ihτ} )} ⩽ Cτ∥a_{ij} ∥_{C(Ω̄ )}∥u∥_{H^{1} ((0,T),H^{2} (Ω))} ,
(5.9)  ∥χ_{i}∥_{L^{2} (Q_{ihτ} )} ⩽ Cτ^{α} ∥u∥_{C_{+}^{α} ([0,T],H^{2} (Ω))} ,
(5.10)  ∥μ_{i}∥_{L^{2} (Q_{ihτ} )} ⩽ Cτ^{2α} ∥u∥_{C_{+}^{α} ([0,T],H^{3} (Ω))} .

The result follows from (5.4)–(5.10). □

In the case when Lu = −Δu we shall derive error bound in the norm H_{+}^{2,α} (Q_{hτ}). Then, for f ∈ C(Q̄), equation (4.1) reduces to

(5.11)  ((I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} v)^{m} + Av^{m−1} = f^{m}, x ∈ ω_{h},

where denoted Av = (A_{1} + A_{2})v = −v_{x_{1} x̄_{1}} − v_{x_{2} x̄_{2}} .

Let u be the solution of initial-boundary-value problem (3.1)–(3.3) with Lu = −Δu, and v the solution of difference problem (5.11), (4.2), (4.3). The error z = u − v satisfies finite the difference scheme

(5.12)  ((I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} z)^{m} + Az^{m−1} = ψ̄^{m} , x ∈ ω_{h},

and homogeneous boundary and initial conditions (5.2)–(5.3). Here

ψ̄^{m} = (I + θτ^{α} A_{1})(I + θτ^{α} A_{2})D_{t,0+,τ}^{α} u^{m} − Δu^{m−1} − f^{m}

= ξ̄^{m} + μ̄^{m} + ∑_{i=1}^{2} (η̄_{i}^{m} + ζ̄_{i}^{m} + χ̄_{i}^{m} ),

where denoted

ξ̄ = D_{t,0+,τ}^{α} u − (D_{t,0+}^{α} u) = ξ_{1}, η_{i} = \frac{∂^{2}u}{∂x_{i}^{2}} − u_{x_{i} x̄_{i}} , ζ_{i} = τ u_{x_{i} x̄_{i} t̄},

χ̄_{i} = −θτ^{α} D_{t,0+,τ}^{α} u_{x_{i} x̄_{i}} = χ_{i,x_{i}} , μ̄ = θ^{2}τ^{2α} D_{t,0+,τ}^{α} u_{x̄_{1} x̄_{1} x_{2} x̄_{2}} = μ_{1,x_{1}} + μ_{2,x_{2}} .

Theorem 5.3. Let the assumptions of Theorem 4.2 hold and let the solution u of initial-boundary value problem (3.1)–(3.3) with Lu = −Δu belong to the space C^{2}([0, T], C(Ω̄)) ∩ C_{+}^{α} ([0, T], H^{4}(Ω)) ∩ H^{1}((0, T), H^{2}(Ω)). Then the solution v of finite difference scheme (5.11), (4.2), (4.3) converges to u and the following convergence rate estimate holds:

∥u − v∥_{H_{+}^{2,α} (Q_{hτ} )} = O(h^{2} + τ^{α}).

Proof. Note that under our assumptions the function f(x, t) is continuous, so (5.11) is well defined.

Applying a priori estimate (4.7) to (5.12) we obtain:

(5.13) ∥u − v∥_{H_{+}^{2,α} (Q_{hτ} )} = ∥z∥_{H_{+}^{2,α} (Q_{hτ} )} ⩽ C∥ψ̄∥_{L^{2} (Q_{hτ} )}

⩽ C[∥ξ̄∥_{L^{2} (Q_{hτ} )} + ∥μ̄∥_{L^{2} (Q_{hτ} )} + ∑_{i=1}^{2} (∥η̄_{i}∥_{L^{2} (Q_{hτ} )} + ∥ζ̄_{i}∥_{L^{2} (Q_{hτ} )} + ∥χ̄_{i} ∥_{L^{2} (Q_{hτ} )} )].

From Lemma 4.1 immediately it follows:

(5.14)  ∥ξ̄∥_{L^{2} (Q_{hτ} )} ⩽ Cτ^{2−α} ∥u∥_{C^{2} ([0,T],C(Ω̄ ))} .

The term η̄_{i} can be estimated using the Bramble–Hilbert lemma:

(5.15)  ∥η̄_{i}∥_{L^{2} (Q_{hτ} )} ⩽ Ch^{2}∥u∥_{C([0,T],H^{4} (Ω))} .

The terms ζ̄_{i}, χ̄_{i} and μ̄ can be estimated directly:

(5.16)  ∥ζ̄_{i}∥_{L^{2} (Q_{hτ} )} ⩽ Cτ∥u∥_{H^{1} ((0,T),H^{2} (Ω))} ,
(5.17)  ∥χ̄_{i} ∥_{L^{2} (Q_{hτ} )} ⩽ Cτ^{α} ∥u∥_{C_{+}^{α} ([0,T],H^{2} (Ω))} ,
(5.18)  ∥μ̄∥_{L^{2} (Q_{hτ} )} ⩽ Cτ^{2α} ∥u∥_{C_{+}^{α} ([0,T],H^{4} (Ω))} .

The result follows from (5.13)–(5.18). □

References

1. R. A. Adams, J. J. F. Fournier, Sobolev Spaces, 2nd ed., Pure and Applied Mathematics Series 140, Elsevier/Academic Press, Amsterdam, 2003.

2. J. H. Bramble, S. R. Hilbert, Estimation of linear functionals on Sobolev spaces with application to Fourier transform and spline interpolation, SIAM J. Numer. Anal. 7 (1970), 112–124.

3. M. Chen, W. Deng, A second-order numerical method for two-dimensional two-sided space fractional convection diffusion equation, arXiv:1304.3788v1 [math.NA], 13 Apr 2013, 1–25.

4. A. Delić, B. S. Jovanović, Numerical approximation of an interface problem for fractional in time diffusion equation, Appl. Math. Comput. 229 (2014), 467–479.

5. V. J. Ervin, J. P. Roop, Variational formulation for the stationary fractional advection dispersion equation, Numer. Methods Partial Differ. Equations 23 (2006), 558–576.

6. S. Hodžić, Factorized difference scheme for 2D fractional in time diffusion equation, Appl. Anal. Discrete Math. 9 (2015), 199–208.

7. S. Hodžić, B. S. Jovanović, Additive difference scheme for two-dimensional fractional in time diffusion equation, Filomat, to appear.

8. B. S. Jovanović, E. Süli, Analysis of Finite Difference Schemes, Springer Ser. Comput. Math. 46, Springer, 2013.

9. A. Kilbas, H. Srivastava, J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier Science and Technology, Boston, 2006.

10. O. A. Ladyzhenskaya, N. N. Ural’tseva, Linear and Quasilinear Elliptic Equations, Nauka, Moscow, 1964 (in Russian); English edition: Academic Press, New York, 1968.

11. L. Li, D. Xu, Alternating direction implicit-Euler method for the two-dimensional fractional evolution equation, J. Comput. Phys. 236 (2013), 157Ű-168.

12. , A space-time spectral method for the time fractional diffusion equation, SIAM J. Numer. Anal. 47 (2009), 2108–2131.

13. J. L. Lions, E. Magenes, Non Homogeneous Boundary Value Problems and Applications, Springer-Verlag, Berlin and New York, 1972.

14. V. G. Maz’ya, T. O. Shaposhnikova, Theory of Multipliers in Spaces of Differentiable Functions, Monographs and Studies in Mathematics 23, Pitman, Boston, Mass., 1985.

15. K. B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.

16. I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999.

17. A. A. Samarskiĭ, Theory of Difference Schemes, Nauka, Moscow, 1989 (in Russian); English edition: Pure Appl. Math. 240, Marcel Dekker, 2001.

18. Z. Z. Sun, X. N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006), 193–209.

19. P. Zhuang, F. Liu, Finite difference approximation for two-dimensional time fractional diffusion equation, J. Algorithms Comput. Technol. 1 (2007), 1–15.

University of Belgrade (Received 14 05 2015) Faculty of Mathematics Belgrade Serbia

 adelic@matf.bg.ac.rs
 sandra@matf.bg.ac.rs
 bosko@matf.bg.ac.rs