# A Mixed Finite Element Formulation for the Conservative Fractional Diffusion Equations

Authors: Suxiang Yang; Huanzhen Chen
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

We consider a boundary-value problem of one-side conservative elliptic equation involving Riemann-Liouville fractional integral. The appearance of the singular term in the solution leads to lower regularity of the solution of the equation, so to the lower order convergence rate for the numerical solution. In this paper, by the dividing of equation, we drop the lower regularity term in the solution successfully and get a new fractional elliptic equation which has full regularity. We present a theoretical framework of mixed finite element approximation to the new fractional elliptic equation and derive the error estimates for unknown function, its derivative, and fractional-order flux. Some numerical results are illustrated to confirm the optimal error estimates.


---

## Body

## 1. Introduction

We consider the following one-side conservative diffusion problems of order2
-
β: (1a)
-
D
K
I
x
β
0
D
u
x
=
f
x
,
x
∈
Ω
≔
0,1
,
(1b)
u
0
=
0
,
u
1
=
0
,where D is the first-order derivative operator, I
x
β
0, defined by (2), refers to the left Riemann-Liouville fractional integral of order β
∈
[
0,1
), K is a positive constant diffusive coefficient, and f
(
x
)
∈
L
2
(
Ω
) is the source or sink term. Equations (1a) and (1b) reduce to the classical second-order diffusion equation when β
=
0.Equations (1a) and (1b), delicately derived from random walk model [1], fractional conservation of mass [2], and the impact of boundary [3], describe subsurface fluid flow and solute transport processes taking place in a multiscale heterogeneous medium/aquifer (see [4] and the references cited therein). These processes exhibit anomalous diffusion or non-Fickian diffusion that cannot be modeled adequately and properly by classical second-order diffusion equation and thus have contracted considerable attention in practical applications.Since the analytic solutions to general fractional partial differential equations are rarely available, one has to turn to numerical methods in general. In the last decade, a number of numerical methods, for example, finite difference method [5], finite volume method [6], spectral method [7], fast difference method [8], finite element method [9–13], and Petrov-Galerkin formulation [14–18], have been developed consecutively for space-fractional partial differential equations.Considering the conservation property of the diffusion processes and the application of the fractional-order fluxp
=
-
K
I
x
β
0
D
u in engineering, an ideal numerical procedure should recognize both the unknown function and its flux and obey the conservation of mass to reflect the physical character of the diffusion model ((1a) and (1b)). Hence, one should involve the state variable u and the fractional-order flux p
=
-
K
I
x
β
0
D
u to form a saddle-point formulation then design a locally conservative numerical procedure. For this purpose, by introducing two intermediate variables q
=
D
u and the flux p
=
-
K
I
x
β
0
D
u, Chen and Wang [19] proposed a saddle-point framework and developed a locally conservative expanded mixed finite element method recently. Although a rigorous numerical analysis theory was established in [19], the convergence rates for the unknown u and the flux p are not optimal with respect to regularity requirement for sufficiently smooth solution u. In particular, for the solution u only in H
1
-
β
+
γ
(
Ω
) with γ
∈
[
β
/
2
,
1
/
2
), the convergence rates are heavily destroyed to γ
-
β
/
2, due to the appearance of the singular term x
1
-
β.The aim of this paper is to express the solution of (1a) and (1b) by a full-regular solution satisfying a general fractional diffusion equation with f-dependent right-hand side term and a x
1
-
β-type solution satisfying an analytic-solved fractional equation. Then, the expanded mixed finite element method proposed in [19] can be applied to solve the f-dependent fractional equation and obtain good error estimates whatever the regularity of the solution u is.The rest of the paper is organized as follows. In Section2, we revisit the fractional differential operators and their properties. In Section 3, we split the fractional diffusion equations (1a) and (1b) into two equations: one is an f-dependent fractional diffusion equation that permits a fully regular solution, and the other is an analytic-solved fractional equation from which the x
1
-
β-type analytic solution is easily solved. By doing so, the solution u of (1a) and (1b) is expressed as the sum of the f-dependent solution and the x
1
-
β-type solution. We are devoted to discretizing the f-dependent fractional equation by the expanded mixed finite element procedure in Section 4 and deriving the related error estimates in Section 5. Obviously, the error estimates for the original fractional diffusion equations (1a) and (1b) obey those for the f-dependent fractional equation since the x
1
-
β-type solution is analytically solved. In Section 6, numerical experiments are performed to confirm our theoretical findings.


## 2. Fractional Operators and Fractional Sobolev Spaces

In this section, we will recall the definitions and properties of Riemann-Liouville fractional differential operators and fractional Sobolev spaces.Definition 1 (see [20–22]).
Letα
>
0. The left and right Riemann-Liouville fractional integrals of order α are defined by (2)
I
x
α
0
u
=
1
Γ
α
∫
0
x
x
-
t
α
-
1
u
t
d
t
,
(3)
I
1
α
x
u
=
1
Γ
α
∫
x
1
t
-
x
α
-
1
u
t
d
t
,where Γ
(
·
) is Gamma function.Definition 2 (see [20–22]).
Letα
>
0 with n
-
1
<
α
≤
n. The left and right Riemann-Liouville fractional derivatives of order α are (4)
D
x
α
0
u
=
d
n
d
x
n
I
x
n
-
α
0
u
,
D
1
α
x
u
=
-
1
n
d
n
d
x
n
I
1
n
-
α
x
u
.We note that fractional integral operator [20–22] satisfies semigroup property; that is, for α
,
β
>
0, (5)
I
x
α
0
I
x
β
0
u
=
I
x
α
+
β
0
u
,
I
1
α
x
I
1
β
x
u
=
I
1
α
+
β
x
u
,
∀
u
∈
L
2
Ω
,and adjoint property (6)
I
x
α
0
u
,
v
=
u
,
I
1
α
x
v
,
∀
u
,
v
∈
L
2
Ω
.We also observe that the left (right) fractional derivative operator is a left inverse of the left (right) fractional integral operator [20–22]: (7)
D
x
α
0
I
x
α
0
u
=
u
,
D
1
α
x
I
1
α
x
u
=
u
,
∀
u
∈
L
2
Ω
.We introduce the fractional derivative spaces and present their well-established equivalence to fractional-order Sobolev spaces [9, 19, 23].Definition 3 (see [9]).
Forμ
>
0, define the norm (8)
u
J
L
μ
R
≔
u
2
+
u
J
L
μ
R
2
1
/
2
,with seminorm (9)
u
J
L
μ
R
≔
D
x
μ
0
u
.Let J
L
μ
(
R
) be the left fractional derivative space defined as the closures of C
0
∞
(
R
) with respect to norm ·
J
L
μ
(
R
).Definition 4 (see [23]).
Forμ
>
0, define the norm (10)
u
J
L
-
μ
R
≔
I
x
μ
0
u
,and let J
L
-
μ
(
R
) denote the closures of C
0
∞
(
R
) with respect to norm ·
J
L
-
μ
(
R
).The right fractional derivative spacesJ
R
μ
(
R
) and J
R
-
μ
(
R
) and their respective norms are defined analogously.Definition 5 (see [9]).
Forμ
>
0, define the norm (11)
u
H
μ
R
≔
u
2
+
u
H
μ
R
2
1
/
2
,with seminorm (12)
u
H
μ
R
≔
ω
μ
u
^
ω
,where u
^
(
ω
) denotes the Fourier transform of u
(13)
u
^
ω
=
∫
R
e
-
i
ω
x
u
x
d
x
.Let H
μ
(
R
) be the left fractional derivative space defined as the closures of C
0
∞
(
R
) with respect to norm ·
H
μ
(
R
).Definition 6 (see [23]).
Forμ
>
0, define the norm (14)
u
H
-
μ
R
≔
ω
-
μ
u
^
ω
,and let H
-
μ
(
R
) denote the closures of C
0
∞
(
R
) with respect to norm ·
H
-
μ
(
R
).We then reiterate the equivalence theories established in [9, 23] for the fractional derivative spaces.Lemma 7 (see [9]).
Letμ
>
0 and μ
≠
n
-
1
/
2, n
=
1,2
,
…. Then, J
L
μ
(
R
), J
R
μ
(
R
), and the fractional-order Sobolev spaces H
μ
(
R
) are equal, with equivalent seminorms and norms.Lemma 8 (see [23]).
Letμ
>
0. Then, J
L
-
μ
(
R
), J
R
-
μ
(
R
), and the negative fractional-order Sobolev spaces H
-
μ
(
R
), the dual space of H
μ
(
R
), are equal with equivalent norms.We now restrict the fractional derivative spaces toΩ.Definition 9 (see [9]).
Letμ
>
0. Define the spaces J
L
,
0
μ
(
Ω
), J
R
,
0
μ
(
Ω
), and H
0
μ
(
Ω
) as the closures of C
0
∞
(
Ω
) under their respective norms.Definition 10 (see [23]).
Letμ
>
0. Define the spaces J
L
-
μ
(
Ω
), J
R
-
μ
(
Ω
), and H
-
μ
(
Ω
) as the closures of C
0
∞
(
Ω
) under their respective norms.Lemma 11 (see [9]).
Letμ
>
0 and μ
≠
n
-
1
/
2, n
=
1,2
,
…. Then, J
L
,
0
μ
(
Ω
), J
R
,
0
μ
(
Ω
), and the fractional-order Sobolev spaces H
0
μ
(
Ω
) are equal, with equivalent seminorms and norms.Lemma 12 (see [19]).
Letμ
>
0. Then, J
L
-
μ
(
Ω
), J
R
-
μ
(
Ω
), and the negative fractional-order Sobolev spaces H
-
μ
(
Ω
), the dual space of H
0
μ
(
Ω
), are equal with equivalent norms.For simplicity of presentation, we use·
μ and ·
μ to stand for the norms and seminorms of fractional-order Sobolev spaces H
μ
(
Ω
) and ·
-
μ to stand for the norms of negative fractional-order Sobolev spaces H
-
μ
(
Ω
), respectively.


## 3. Decomposition of (1a) and (1b)

In this section, we shall split the fractional diffusion equations (1a) and (1b) into an f-dependent fractional diffusion equation that permits a fully regular solution and an analytic-solved fractional equation from which the x
1
-
β-type analytic solution is easily solved.Theorem 13.
Letu be the solution to the fractional diffusion problem ((1a) and (1b)) with right term f
(
x
)
∈
L
2
(
Ω
). Then,(1)
u can be expressed as (15)
u
=
u
r
+
u
s
,where u
r and u
s are the solutions to the following fractional problems (16a)
-
D
K
I
x
β
0
D
u
r
x
=
f
r
,
x
∈
Ω
,
(16b)
u
r
0
=
0
,
u
r
1
=
0
,
(17a)
-
D
K
I
x
β
0
D
u
s
x
=
2
K
M
0
Γ
1
+
β
x
β
,
x
∈
Ω
,
(17b)
u
s
0
=
0
,
u
s
1
=
0
;here f
r
=
f
(
x
)
-
2
K
M
0
/
Γ
(
1
+
β
)
x
β with M
0
=
I
x
0
D
x
β
0
I
x
0
(
f
/
K
)
(
1
);
(2)
letp
≔
-
K
I
x
β
0
D
u and p
r
≔
-
K
I
x
β
0
D
u
r be fractional-order flux. Then, u
∈
H
0
1
-
β
+
γ
(
Ω
) with β
/
2
≤
γ
<
1
/
2, p
∈
H
1
(
Ω
), u
r
∈
H
0
2
-
β
(
Ω
), p
r
∈
H
1
(
Ω
); further, (18)
u
1
-
β
+
γ
+
p
1
≤
C
f
,
u
r
2
-
β
+
p
r
1
≤
C
f
.Proof.
The solution and the regularity of the fractional diffusion equations (1a) and (1b) are given in [19], and the proof is provided only for completeness. Introducing two intermediate variables q
≔
D
u and p
≔
-
K
I
x
β
0
q, ((1a) and (1b)) can be rewritten as (19a)
D
u
=
q
,
(19b)
I
x
β
0
q
=
-
p
K
,
(19c)
D
p
=
f
.Solving (19c) to obtain (20)
p
=
p
0
+
I
x
0
fand noting I
x
0
f
∈
H
1
(
Ω
) for f
∈
L
2
(
Ω
), we conclude that p
∈
H
1
(
Ω
) and p
1
≤
C
f.
Solving (19a) and applying D
x
β
0 to both sides of (19b), we deduce (21)
u
=
u
0
+
I
x
0
q
=
u
0
-
I
x
0
D
x
β
0
p
K
.Substituting (20) into the above equation and noting the boundary condition u
(
0
)
=
0, we get (22)
u
=
-
p
0
K
Γ
2
-
β
x
1
-
β
-
I
x
0
D
x
β
0
I
x
0
f
K
.Enforcing the boundary condition u
(
1
)
=
0 to the above equation, we obtain -
p
(
0
)
/
K
Γ
(
2
-
β
)
=
I
x
0
D
x
β
0
I
x
0
(
f
/
K
)
(
1
)
≔
M
0, which leads to the solution of (1a) and (1b):(23)
u
=
M
0
x
1
-
β
-
I
x
0
D
x
β
0
I
x
0
f
K
.
Next, we are to solve (16a) and (16b) and (17a) and (17b). Substituting f
(
x
)
-
2
K
M
0
/
Γ
(
1
+
β
)
x
β for f
(
x
) in (22) and noticing the boundary condition u
r
(
1
)
=
0, we get the solution of (16a) and (16b): (24)
u
r
x
=
-
I
x
0
D
x
β
0
I
x
0
f
K
+
M
0
x
2
.Similarly, we get the the solution of (17a) and (17b): (25)
u
s
x
=
M
0
x
1
-
β
-
x
2
.It follows from (23)–(25) that the solution to the fractional problem ((1a) and (1b)) can be split into the solutions of (16a) and (16b) and (17a) and (17b).
Noting thatD
x
β
0 is a bounded operator from fractional Sobolev space H
β
(
Ω
) to L
2
(
Ω
) [3], we get the first term -
I
x
0
D
x
β
0
I
x
0
(
f
/
K
) of solution u in (23) is in H
2
-
β
(
Ω
) with f
∈
L
2
(
Ω
). The second term M
0
x
1
-
β
∈
H
1
-
β
+
γ
(
Ω
) with β
/
2
≤
γ
<
1
/
2 [14]. Combining the regularities of two parts and the boundary conditions, we derive u
∈
H
0
1
-
β
+
γ
(
Ω
) and u
1
-
β
+
γ
≤
C
f. Clearly, the regularity of the second term M
0
x
2
∈
H
2
+
γ
(
Ω
) is no less than the regularity of the first term in (24). This implies u
r
∈
H
0
2
-
β
(
Ω
), p
r
∈
H
1
(
Ω
), and u
r
2
-
β
+
p
r
1
≤
C
f.Remark 14.
The decomposition of (1a) and (1b) is not unique in Theorem 13; that is, we can choose any other function v, provided that the regularity of v is no less than the regularity of the first term in (24).Remark 15.
From Theorem13, the regularity of original problem (1a) and (1b) is only H
0
1
-
β
+
γ
(
Ω
), due to the presence of the singular term x
1
-
β, while the regularity to (18) is H
0
2
-
β
(
Ω
) with f
(
x
)
∈
L
2
(
Ω
). In particular, for the right term f
(
x
)
∈
H
l
(
Ω
), l
≥
0, the regularity can be improved to H
0
l
+
2
-
β
(
Ω
) by choosing a proper decomposition. For example, we can replace x
2 by x
m with m
≥
l
+
2
-
β in the decomposition of the solution when f
(
x
)
∈
H
l
(
Ω
). This is the advantage of the splitting.Remark 16.
Equations (17a) and (17b) is analytic-solvable. Once the regular part (
p
r
,
q
r
,
u
r
) is determined, the solution (
p
,
q
,
u
) of the source problem ((19a), (19b), and (19c)) (or (1a) and (1b)) has the representation (26)
p
=
p
r
-
K
M
0
Γ
2
-
β
-
2
Γ
2
+
β
x
1
+
β
,
q
=
q
r
+
M
0
1
-
β
x
-
β
-
2
x
,
u
=
u
r
+
M
0
x
1
-
β
-
x
2
.So we focus on f-dependent fractional diffusion equations (16a) and (16b) from now on.


## 4. Saddle-Point Formulation and Mixed Finite Element Method

In this section, we shall formulate a saddle-point variational formulation and establish a mixed finite element formulation to thef-dependent fractional diffusion equations (16a) and (16b).Introduce the fractional-order fluxp
r
=
-
K
I
x
β
0
q
r with q
r
=
D
u
r. Then, we can rewrite (16a) and (16b) into the following mixed form: (27a)
D
u
r
=
q
r
,
(27b)
I
x
β
0
q
r
=
-
p
r
K
,
(27c)
D
p
r
=
f
r
.The mixed variational formulation corresponding to (27a), (27b), and (27c) is to find (
p
r
,
q
r
,
u
r
)
∈
V
×
H
×
W such that (28a)
q
r
,
v
+
u
r
,
D
v
=
0
,
v
∈
V
,
(28b)
p
r
,
σ
+
K
I
x
β
/
2
0
q
r
,
I
1
β
/
2
x
σ
=
0
,
σ
∈
H
,
(28c)
D
p
r
,
w
=
f
r
,
w
,
w
∈
W
,where V
=
H
1
(
Ω
), H
=
H
-
β
/
2
(
Ω
), and W
=
L
2
(
Ω
).Lemma 17 (see [19]).
There exists a unique solution(
p
r
,
q
r
,
u
r
)
∈
V
×
H
×
W to the problem (28a), (28b), and (28c).Next, we shall establish a mixed finite element procedure based on the saddle-point variational formulation (28a), (28b), and (28c).LetI
h be a uniform division of Ω by I
i
=
[
x
i
-
1
,
x
i
], i
=
1,2
,
…
,
M with x
0
=
0, x
M
=
1, and h
=
1
/
M. Let V
h
×
H
h
×
W
h
⊂
V
×
H
×
W denote Raviart-Thomas space [24–27]: (29)
V
h
=
v
h
∈
V
∩
C
0
Ω
¯
:
v
h
I
i
∈
P
k
I
i
,
i
=
1,2
,
…
,
M
,
k
≥
1
,
H
h
=
σ
h
∈
H
:
σ
h
I
i
∈
P
k
I
i
,
i
=
1,2
,
…
,
M
,
k
≥
1
,
W
h
=
w
h
∈
W
:
w
h
I
i
∈
P
k
-
1
I
i
,
i
=
1,2
,
…
,
M
,
k
≥
1
,where P
k
(
I
i
) is the restriction of the set of all polynomials of a degree not bigger than k on I
i. Then, the mixed finite element formulation is to find (
p
h
r
,
q
h
r
,
u
h
r
)
∈
V
h
×
H
h
×
W
h such that (30a)
q
h
r
,
v
h
+
u
h
r
,
D
v
h
=
0
,
v
h
∈
V
h
,
(30b)
p
h
r
,
σ
h
+
K
I
x
β
/
2
0
q
h
r
,
I
1
β
/
2
x
σ
h
=
0
,
σ
h
∈
H
h
,
(30c)
D
p
h
r
,
w
h
=
f
r
,
w
h
,
w
h
∈
W
h
.Last, we construct an approximate solution (
p
h
,
q
h
,
u
h
) to (19a), (19b), and (19c) (or (1a) and (1b)) by (31)
p
h
=
p
h
r
-
K
M
0
Γ
2
-
β
-
2
Γ
2
+
β
x
1
+
β
,
q
h
=
q
h
r
+
M
0
1
-
β
x
-
β
-
2
x
,
u
h
=
u
h
r
+
M
0
x
1
-
β
-
x
2
.


## 5. Convergence Analysis

In this section, we are devoted to the convergence analysis. In error analysis below, we shall make use of three projection operators. The first operator isΠ
h
:
V
↦
V
h defined by (32)
D
v
-
Π
h
v
,
w
h
=
0
,
∀
w
h
∈
W
h
.The other two operators are the L
2-projections P
h
:
W
↦
W
h and R
h
:
H
↦
H
h: (33)
u
-
P
h
u
,
w
h
=
0
,
∀
w
h
∈
W
h
,
q
-
R
h
q
,
v
h
=
0
,
∀
v
h
∈
H
h
.Lemmas 18–20 hold in the integer-order Sobolev spaces (if n
=
1 see [19] and if n
=
2,3 see [24–27]). We can prove that they also hold in the fractional-order spaces by applying standard real interpolation methods [26].Lemma 18.
Letl and s be any real numbers. Assume that u
∈
H
s
(
Ω
). Then, (34)
u
-
P
h
u
-
l
≤
C
h
s
+
l
u
s
,
0
≤
l
,
s
≤
k
.Lemma 19.
Letl and s be any real numbers. Assume that q
∈
H
s
(
Ω
). Then, (35)
q
-
R
h
q
-
l
≤
C
h
s
+
l
q
s
,
0
≤
l
,
s
≤
k
+
1
.Lemma 20.
Assume thatv
∈
H
s
(
Ω
); then, (36)
v
-
Π
h
v
θ
≤
C
h
s
-
θ
v
s
,
0
≤
θ
≤
1
.Now we proceed to the convergence analysis. Letξ
=
Π
h
p
r
-
p
h
r, η
=
R
h
q
r
-
q
h
r, and ζ
=
P
h
u
r
-
u
h
r; then, from (28a), (28b), and (28c), (30a), (30b), and (30c), (32), and (33), we can get the error equations (37a)
η
,
v
h
+
ζ
,
D
v
h
=
0
,
v
h
∈
V
h
,
(37b)
p
r
-
p
h
r
,
σ
h
+
K
I
x
β
/
2
0
q
r
-
q
h
r
,
I
1
β
/
2
x
σ
h
=
0
,
σ
h
∈
H
h
,
(37c)
D
ξ
,
w
h
=
0
,
w
h
∈
W
h
.Lemma 21.
Suppose thatξ satisfies (37c). Then, the following holds: (38)
D
p
r
-
p
h
r
=
D
p
r
-
Π
h
p
r
.Proof.
Setw
h
=
D
ξ
∈
W
h in (37c) to derive that D
ξ
=
0, which implies ξ
=
constant. Then, we obtain the desired result.Lemma 22.
There exists a constantC
>
0, such that (39)
q
r
-
q
h
r
H
≤
C
q
r
-
R
h
q
r
H
+
C
h
1
-
β
/
2
D
p
r
-
Π
h
p
r
.Proof.
Setv
h
=
ξ in (37a). Then, noting D
ξ
=
0, we have (40)
ξ
,
η
=
0
.Set σ
h
=
η in (37b); then, (41)
C
0
η
H
2
≤
K
cos
⁡
π
β
2
I
x
β
/
2
-
∞
η
~
=
K
I
x
β
/
2
0
η
,
I
1
β
/
2
x
η
=
-
K
I
x
β
/
2
0
q
r
-
R
h
q
r
,
I
1
β
/
2
x
η
-
p
r
-
Π
h
p
r
,
η
,where (
I
x
β
/
2
0
η
,
I
1
β
/
2
x
η
)
=
(
I
x
β
/
2
-
∞
η
~
,
I
∞
β
/
2
x
η
~
)
=
cos
⁡
(
π
β
/
2
)
I
x
β
/
2
-
∞
η
~, η
~ being the extension of η by zeros to R [23, Lemma  2.6], is used and C
0
=
K
cos
⁡
(
π
β
/
2
). Now, we estimate the right-hand side of (41) term by term. For the first term, we get (42)
-
K
I
x
β
/
2
0
q
r
-
R
h
q
r
,
I
1
β
/
2
x
η
≤
C
q
r
-
R
h
q
r
H
η
H
.For the second term, applying Green’s formula, we have (43)
p
r
-
Π
h
p
r
,
η
=
p
r
-
Π
h
p
r
,
D
I
x
0
η
=
-
D
p
r
-
Π
h
p
r
,
I
x
0
η
,where I
x
0
η
(
0
)
=
I
x
0
η
(
1
)
=
0 is used. In fact, taking v
h
=
1 in (37a), we get (
η
,
1
)
=
0; that is, I
x
0
η
(
1
)
=
0. Introduce the L
2-projection operator P
h
0. Then, the approximation property I
x
0
η
-
P
h
0
(
I
x
0
η
)
≤
C
h
1
-
β
/
2
I
x
0
η
1
-
β
/
2 holds. By the property of the L
2-projection operator P
h
0, (43) can be rewritten as (44)
p
r
-
Π
h
p
r
,
η
=
-
D
p
r
-
Π
h
p
r
,
I
x
0
η
-
P
h
0
I
x
0
η
≤
D
p
r
-
Π
h
p
r
I
x
0
η
-
P
h
0
I
x
0
η
≤
C
h
1
-
β
/
2
D
p
r
-
Π
h
p
r
I
x
0
η
1
-
β
/
2
≤
C
h
1
-
β
/
2
D
p
r
-
Π
h
p
r
η
H
,where I
x
0
η
1
-
β
/
2
≤
C
D
x
1
-
β
/
2
0
I
x
0
η
=
C
I
x
β
/
2
0
η
≤
C
η
H is used. Substituting (42) and (44) into (41) yields (45)
η
H
≤
C
h
1
-
β
/
2
D
p
r
-
Π
h
p
r
+
C
q
r
-
R
h
q
r
H
.In conjunction with triangle inequality, we obtain the desired estimate.Lemma 23.
There exists a constantC
>
0, such that (46)
p
r
-
p
h
r
≤
C
Π
h
p
r
-
p
r
+
q
r
-
q
h
r
H
.Proof.
Takingσ
h
=
ξ in error equations (37a), (37b), and (37c), we have (47)
ξ
2
=
Π
h
p
r
-
p
r
,
ξ
+
p
r
-
p
h
r
,
ξ
=
Π
h
p
r
-
p
r
,
ξ
-
K
I
x
β
/
2
0
q
r
-
q
h
r
,
I
1
β
/
2
x
ξ
≤
Π
h
p
r
-
p
r
ξ
+
C
q
r
-
q
h
r
H
ξ
H
≤
Π
h
p
r
-
p
r
ξ
+
C
q
r
-
q
h
r
H
ξ
.Then, (48)
ξ
≤
Π
h
p
r
-
p
r
+
C
q
r
-
q
h
r
H
.Together with triangle inequality, we finish the proof of Lemma 23.Now let us apply the duality argument of Douglas and Roberts to obtain an error estimate in theL
2-norm.Lemma 24.
There exists a constantC
>
0, such that (49)
u
r
-
u
h
r
≤
u
r
-
P
h
u
r
+
C
h
γ
-
β
/
2
q
r
-
q
h
r
H
+
C
h
γ
-
β
p
r
-
p
h
r
+
C
h
1
-
β
+
γ
D
p
r
-
p
h
r
.Proof.
Letg
∈
L
2
(
Ω
) and w
∈
H
0
1
-
β
+
γ
(
Ω
), β
/
2
≤
γ
<
1
/
2, such that (50a)
-
D
K
I
1
β
x
D
w
x
=
g
x
,
x
∈
Ω
=
0,1
,
(50b)
w
0
=
0
,
w
1
=
0
.Then, the regularity estimate w
1
-
β
+
γ
≤
C
g holds. Taking v
h
=
Π
h
(
K
I
1
β
x
D
w
) in (37a), we have (51)
ζ
,
g
=
ζ
,
-
D
K
I
1
β
x
D
w
=
ζ
,
-
D
Π
h
K
I
1
β
x
D
w
=
q
r
-
q
h
r
,
Π
h
K
I
1
β
x
D
w
=
q
r
-
q
h
r
,
K
I
1
β
x
D
w
+
q
r
-
q
h
r
,
Π
h
K
I
1
β
x
D
w
-
K
I
1
β
x
D
w
≔
G
1
+
G
2
.First, we are to estimate the first term G
1 of (51). Using adjoint property (6), we have (52)
G
1
=
q
r
-
q
h
r
,
K
I
1
β
x
D
w
=
K
I
x
β
0
q
r
-
q
h
r
,
D
w
=
K
I
x
β
0
q
r
-
q
h
r
,
D
w
-
R
h
D
w
+
K
I
x
β
0
q
r
-
q
h
r
,
R
h
D
w
≡
G
11
+
G
12
.We use Lemma 19 to estimate G
11: (53)
G
11
≤
q
r
-
q
h
r
H
D
w
-
R
h
D
w
H
≤
C
h
γ
-
β
/
2
D
w
γ
-
β
q
r
-
q
h
r
H
≤
C
h
γ
-
β
/
2
q
r
-
q
h
r
H
w
1
-
β
+
γ
.The estimate of G
12 is the result of (37c), Lemmas 18 and 19: (54)
G
12
=
-
p
r
-
p
h
r
,
R
h
D
w
=
p
r
-
p
h
r
,
D
w
-
R
h
D
w
+
p
r
-
p
h
r
,
D
w
=
p
r
-
p
h
r
,
D
w
-
R
h
D
w
-
D
p
r
-
p
h
r
,
w
=
p
r
-
p
h
r
,
D
w
-
R
h
D
w
-
D
p
r
-
p
h
r
,
w
-
P
h
w
-
D
p
r
-
p
h
r
,
P
h
w
≤
p
r
-
p
h
r
D
w
-
R
h
D
w
+
D
p
r
-
p
h
r
w
-
P
h
w
≤
C
h
γ
-
β
p
r
-
p
h
r
w
1
-
β
+
γ
+
C
h
1
-
β
+
γ
D
p
r
-
p
h
r
w
1
-
β
+
γ
.Together with above two inequalities, we get the estimate of G
1: (55)
G
1
≤
C
h
γ
-
β
/
2
q
r
-
q
h
r
H
+
h
γ
-
β
p
r
-
p
h
r
+
h
1
-
β
+
γ
D
p
r
-
p
h
r
w
1
-
β
+
γ
.
Thanks to the adjoint property of fractional integral operator, the relationship ofL
2-projection operator P
h and R
-
T operator Π
h, and adjoint problem (50a) and (50b), we get the estimate of G
2: (56)
G
2
=
q
r
-
q
h
r
,
Π
h
K
I
1
β
x
D
w
-
K
I
1
β
x
D
w
=
D
I
x
0
q
r
-
q
h
r
,
Π
h
K
I
1
β
x
D
w
-
K
I
1
β
x
D
w
=
-
I
x
0
q
r
-
q
h
r
,
D
Π
h
K
I
1
β
x
D
w
-
K
I
1
β
x
D
w
=
-
I
x
0
q
r
-
q
h
r
,
P
h
D
K
I
1
β
x
D
w
-
D
K
I
1
β
x
D
w
=
I
x
0
q
r
-
q
h
r
,
P
h
g
-
g
=
-
I
x
β
/
2
0
q
r
-
q
h
r
,
I
1
1
-
β
/
2
x
P
h
g
-
g
≤
q
r
-
q
h
r
H
P
h
g
-
g
-
1
-
β
/
2
≤
C
h
1
-
β
/
2
g
q
r
-
q
h
r
H
,where I
x
0
(
q
r
-
q
h
r
)
(
0
)
=
I
x
0
(
q
r
-
q
h
r
)
(
1
)
=
0 is used. In fact, taking v
h
=
1 in (37a), we derive (
q
r
-
q
h
r
,
1
)
=
0; that is, I
x
0
(
q
r
-
q
h
r
)
(
1
)
=
0. Substituting (55) and (56) into (51), and noticing the regularity of the adjoint problem (50a) and (50b), we obtain (57)
ζ
≤
C
h
γ
-
β
/
2
q
r
-
q
h
r
H
+
h
γ
-
β
p
r
-
p
h
r
+
h
1
-
β
+
γ
D
p
r
-
p
h
r
.By triangle inequality, we get the desired result.By Lemmas21–24, the properties of Raviart-Thomas operator Π
h, L
2 projections P
h, and R
h, we obtain the following.Theorem 25.
Letf
(
x
)
∈
H
l
(
Ω
). Suppose u
r is the solution of (27a), (27b), and (27c) and (
p
h
r
,
q
h
r
,
u
h
r
)
∈
V
h
×
H
h
×
W
h satisfies (30a), (30b), and (30c). Then, we can get the following error estimates: (58)
D
p
r
-
p
h
r
≤
C
h
min
⁡
l
,
k
f
l
,
q
r
-
q
h
r
H
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
,
p
r
-
p
h
r
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
,
u
r
-
u
h
r
≤
C
h
min
⁡
l
+
1
-
3
β
/
2
+
γ
,
k
f
l
.Proof.
Letf
(
x
)
∈
H
l
(
Ω
). By the regularity of problem (16a) and (16b), we know u
r
∈
H
l
+
2
-
β
(
Ω
), q
r
=
D
u
r
∈
H
l
+
1
-
β
(
Ω
), and p
r
∈
H
l
+
1
(
Ω
). Recalling the approximation properties in Lemmas 19 and 20, (59)
q
r
-
R
h
q
r
H
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
+
β
/
2
q
r
l
+
1
-
β
,
D
p
r
-
Π
h
p
r
≤
C
h
min
⁡
l
,
k
p
r
l
+
1
,and combining with Lemmas 21 and 22, we can get the estimates of D
(
p
r
-
p
h
r
) and q
r
-
q
h
r
H: (60)
q
r
-
q
h
r
H
≤
C
q
r
-
R
h
q
r
H
+
C
h
1
-
β
/
2
D
p
r
-
Π
h
p
r
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
+
β
/
2
q
r
l
+
1
-
β
+
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
p
r
l
+
1
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
q
r
l
+
1
-
β
+
p
r
l
+
1
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
.
Analogously, combining the above estimate ofq
r
-
q
h
r
H and the estimate of p
r
-
Π
h
p
r, (61)
p
r
-
Π
h
p
r
≤
C
h
min
⁡
l
+
1
,
k
+
1
p
r
l
+
1
,we get from Lemma 23
(62)
p
r
-
p
h
r
≤
C
Π
h
p
r
-
p
r
+
q
r
-
q
h
r
H
≤
C
h
min
⁡
l
+
1
,
k
+
1
p
r
l
+
1
+
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
q
r
l
+
1
-
β
+
p
r
l
+
1
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
q
r
l
+
1
-
β
+
p
r
l
+
1
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
.The estimate of u
r
-
u
h
r is a result of Lemmas 21–24 and the following estimate: (63)
u
r
-
P
h
u
r
≤
C
h
min
⁡
l
+
2
-
β
,
k
u
r
l
+
2
-
β
.By (26), (31), and Theorem 25, we can derive error estimates for the approximation (
p
h
,
q
h
,
u
h
).Theorem 26.
Letf
(
x
)
∈
H
l
(
Ω
). Suppose (
p
,
q
,
u
) is the solution of (19a), (19b), and (19c) (or (1a) and (1b)) and (
p
h
,
q
h
,
u
h
) satisfies (31). Then, we can get the following error estimates: (64)
D
p
-
p
h
≤
C
h
min
⁡
l
,
k
f
l
,
q
-
q
h
H
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
,
p
-
p
h
≤
C
h
min
⁡
l
+
1
-
β
/
2
,
k
+
1
-
β
/
2
f
l
,
u
-
u
h
≤
C
h
min
⁡
l
+
1
-
3
β
/
2
+
γ
,
k
f
l
.Remark 27.
From Theorem26, the estimates D
(
p
-
p
h
) and q
-
q
h
H are optimal, and the estimates of p
-
p
h and u
-
u
h are suboptimal, even though the regularity of the solution u remains only in H
1
-
β
+
γ
(
Ω
) for γ
∈
[
β
/
2
,
1
/
2
). Hence, in comparison with the standard mixed finite element methods [19], the error rates for q
-
q
h
H, p
-
p
h, and u
-
u
h are only γ
-
β
/
2; in spite of the enough smoothness of the right term f, the new approach does yield a higher-order convergence rate.Remark 28.
The estimates ofp
-
p
h and u
-
u
h are only suboptimal, due to the structure of the error equations (37a), (37b), and (37c) and the low global regularity of the adjoint problem (50a) and (50b), used in Nitsches trick. The error equations (37a), (37b), and (37c) involve u
r
,
p
r, and q
r as a whole, so the estimates for u
r
-
u
h
r and p
r
-
p
h
r heavily depend on q
r
-
q
h
r
H. This may weaken the orders of error estimates.


## 6. Numerical Experiments

In this part, we present some numerical experiments to verify the theoretically proven convergence results. In all the following examples, we letK
=
1. For each example, we consider three different β values, that is, 1
/
4, 1
/
2, and 2
/
3, and present the estimates and the convergence orders of u
-
u
h, p
-
p
h, and q
-
q
h
H
-
β
/
2
(
Ω
), where the number in the bracket under the column order is the theoretical convergence rate derived in Theorem 26.We consider four examples with the analytic solution of different regularity. The analytic solution of Example1 is sufficiently smooth, while the others only belong to H
1
-
β
+
γ
(
Ω
) with γ
∈
[
β
/
2
,
1
/
2
), due to the appearance of singular term x
1
-
β
∈
H
1
-
β
+
γ
(
Ω
). So we split (1a) and (1b) into a f-dependent fractional diffusion equations (16a) and (16b) and an analytic-solved fractional equation (17a) and (17b) by Theorem 13. Spaces V
h
×
H
h
×
W
h
⊂
V
×
H
×
W choose to be (29) with k
=
1 and k
=
2 for each example.Example 1.
Letϵ
∈
[
0,1
/
2
); source term f
(
x
)
=
-
2
x
β
/
Γ
(
β
+
1
)
+
6
x
β
+
1
/
Γ
(
β
+
2
)
∈
H
β
+
ϵ
(
Ω
). Then, the analytic solution u
(
x
)
=
x
2
(
1
-
x
). As u
∈
H
2
+
ϵ
(
Ω
) is very smooth, the result of Theorem 26 predicts almost O
(
h
min
⁡
(
2
-
β
/
2
,
k
)
), O
(
h
min
⁡
3
/
2
+
β
/
2
,
k
+
1
-
β
/
2
), and O
(
h
min
⁡
(
3
/
2
+
β
/
2
,
k
+
1
-
β
/
2
)
) for u
-
u
h
,
p
-
p
h, and q
-
q
h
H
-
β
/
2
(
Ω
), respectively. Tables 1 and 2 include numerical results for k
=
1 and k
=
2, respectively. From Tables 1 and 2, we can see that they are much higher than the theoretical prediction of Theorem 26.Table 1
Numerical result of Example1 on a uniform mesh of size h for k
=
1.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
2.57e − 2
6.58e − 3
1.65e − 3
4.12e − 4
0.99 (1.00)
p
-
p
h
1.90e − 2
1.45e − 3
1.14e − 4
9.36e − 6
1.83 (1.63)
q
-
q
h
H
1.55e − 2
1.35e − 3
1.24e − 4
1.21e − 5
1.72 (1.63)
1/2
u
-
u
h
2.64e − 2
6.60e − 3
1.65e − 3
4.12e − 4
1.00 (1.00)
p
-
p
h
1.72e − 2
1.22e − 3
7.89e − 5
5.01e − 6
1.93 (1.75)
q
-
q
h
H
2.41e − 2
2.13e − 3
1.60e − 4
1.18e − 5
1.77 (1.75)
2/3
u
-
u
h
2.70e − 2
6.61e − 3
1.65e − 3
4.12e − 4
1.01 (1.00)
p
-
p
h
1.72e − 2
1.22e − 3
7.90e − 5
5.01e − 6
1.96 (1.67)
q
-
q
h
H
2.41e − 2
2.13e − 3
1.61e − 4
1.18e − 5
1.83 (1.67)Table 2
Numerical result of Example1 on a uniform mesh of size h for k
=
2.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
7.56e − 3
4.61e − 4
2.88e − 5
1.80e − 6
2.01 (1.88)
p
-
p
h
2.95e − 2
1.52e − 3
6.92e − 5
3.10e − 6
2.20 (1.63)
q
-
q
h
H
4.44e − 2
2.79e − 3
1.52e − 4
8.12e − 6
2.07 (1.63)
1/2
u
-
u
h
7.62e − 3
4.61e − 4
2.88e − 5
1.80e − 6
2.00 (1.75)
p
-
p
h
1.30e − 2
3.89e − 4
1.21e − 5
3.81e − 7
2.51 (1.75)
q
-
q
h
H
3.12e − 2
1.48e − 3
6.70e − 5
3.04e − 6
2.22 (1.75)
2/3
u
-
u
h
7.94e − 3
4.62e − 4
2.88e − 5
1.80e − 6
2.01 (1.67)
p
-
p
h
8.73e − 3
1.59e − 4
3.62e − 6
9.66e − 8
2.74 (1.83)
q
-
q
h
H
2.75e − 2
9.84e − 4
3.86e − 5
1.59e − 6
2.34 (1.83)Example 2.
Let source termf
(
x
)
=
x
1
+
β
∈
H
1
+
β
+
ϵ
(
Ω
) for ϵ
∈
[
0,1
/
2
), and c
β
=
Γ
(
2
+
β
)
/
6. Then, the analytic solution u
(
x
)
=
c
β
(
x
1
-
β
-
x
3
) only belongs to H
1
-
β
+
γ
(
Ω
) with γ
∈
[
β
/
2
,
1
/
2
), despite the smoothness of the right term f
(
x
). So we split (1a) and (1b) into an f-depended fractional diffusion equation with right term f
r
(
x
)
=
x
1
+
β
(
1
-
4
x
/
2
+
β
) and analytic solution u
r
(
x
)
=
c
β
(
x
4
-
x
3
)
∈
H
3
+
ϵ
(
Ω
), and an analytic-solved equation. The numerical results are shown in Tables 3 and 4. The convergence rates are identical with that for Example 1. This is attributed to the good regularity of the analytic solution u
r.Table 3
Numerical result of Example2 on a uniform mesh of size h for k
=
1.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
9.60e − 2
2.42e − 2
6.05e − 3
1.51e − 3
0.98 (1.00)
p
-
p
h
2.06e − 2
1.29e − 3
8.09e − 5
5.06e − 6
2.00 (1.88)
q
-
q
h
H
1.35e − 2
8.15e − 4
4.72e − 5
2.71e − 6
2.05 (1.88)
1/2
u
-
u
h
9.61e − 2
2.42e − 2
6.05e − 3
1.51e − 3
1.00 (1.00)
p
-
p
h
1.95e − 2
1.24e − 3
7.78e − 5
4.87e − 6
1.99 (1.75)
q
-
q
h
H
1.28e − 2
9.43e − 4
6.14e − 5
3.88e − 6
1.95 (1.75)
2/3
u
-
u
h
9.62e − 2
2.42e − 2
6.05e − 3
1.51e − 3
0.99 (1.00)
p
-
p
h
1.89e − 2
1.23e − 3
7.71e − 5
4.82e − 6
1.99 (1.67)
q
-
q
h
H
1.74e − 2
1.36e − 3
8.92e − 5
5.68e − 6
1.93 (1.67)Table 4
Numerical result of Example2 on a uniform mesh of size h for k
=
2.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
1.29e − 2
7.98e − 4
4.98e − 5
3.11e − 6
2.00 (2.00)
p
-
p
h
2.97e − 2
1.52e − 3
6.95e − 5
3.10e − 6
2.20 (2.63)
q
-
q
h
H
4.46e − 2
2.80e − 3
1.52e − 4
8.08e − 6
2.07 (2.63)
1/2
u
-
u
h
1.29e − 2
7.98e − 4
4.98e − 5
3.11e − 6
2.00 (2.00)
p
-
p
h
1.29e − 2
3.86e − 4
1.19e − 5
3.71e − 7
2.51 (2.75)
q
-
q
h
H
3.11e − 2
1.48e − 3
6.64e − 5
2.95e − 6
2.23 (2.75)
2/3
u
-
u
h
1.31e − 2
7.98e − 4
4.98e − 5
3.11e − 6
2.01 (2.00)
p
-
p
h
8.43e − 3
1.49e − 4
3.24e − 6
7.77e − 8
2.79 (2.67)
q
-
q
h
H
2.70e − 2
9.69e − 4
3.76e − 5
1.47e − 6
2.36 (2.67)Example 3.
Let source termf
(
x
)
=
1
∈
H
ϵ
(
Ω
) for ϵ
∈
[
0,1
/
2
), and constants c
β
=
1
/
Γ
(
3
-
β
). The analytic solution u
(
x
)
=
c
β
(
x
1
-
β
-
x
2
-
β
) belongs to H
1
-
β
+
γ
(
Ω
) with γ
∈
[
β
/
2
,
1
/
2
). We split (1a) and (1b) into (16a) and (16b) and (17a) and (17b). We are to solve (16a) and (16b) with right term f
r
(
x
)
=
1
-
2
c
β
/
Γ
(
1
+
β
)
x
β. The analytic solution u
r
(
x
)
=
c
β
(
x
2
-
x
2
-
β
)
∈
H
2
-
β
+
ϵ
(
Ω
) is very smooth. The numerical results are shown in Tables 5 and 6. The same convergence rates are observed, which are much higher than the theoretical prediction of Theorem 26.Table 5
Numerical result of Example3 on a uniform mesh of size h for k
=
1.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
4.95e − 2
1.24e − 2
3.10e − 3
7.76e − 4
1.00 (1.00)
p
-
p
h
1.88e − 6
1.63e − 7
2.50e − 9
3.85e − 11
2.60 (1.38)
q
-
q
h
H
3.57e − 3
5.33e − 4
7.93e − 5
1.18e − 5
1.37 (1.38)
1/2
u
-
u
h
5.74e − 2
1.44e − 2
3.60e − 3
9.00e − 4
1.00 (1.00)
p
-
p
h
1.54e − 4
3.30e − 7
5.60e − 9
8.68e − 11
3.46 (1.25)
q
-
q
h
H
6.70e − 3
1.20e − 3
2.13e − 4
3.76e − 5
1.25 (1.25)
2/3
u
-
u
h
6.25e − 2
1.56e − 2
3.91e − 3
9.78e − 4
1.00 (1.00)
p
-
p
h
6.56e − 4
1.20e − 6
5.67e − 9
1.06e − 10
3.76 (1.17)
q
-
q
h
H
7.67e − 3
1.63e − 3
3.24e − 4
6.44e − 5
1.15 (1.17)Table 6
Numerical result of Example3 on a uniform mesh of size h for k
=
2.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
4.00e − 3
2.58e − 4
1.64e − 5
1.03e − 6
1.99 (1.63)
p
-
p
h
3.88e − 3
5.02e − 4
6.30e − 5
7.88e − 6
1.49 (1.38)
q
-
q
h
H
6.11e − 3
9.41e − 4
1.40e − 4
2.09e − 5
1.37 (1.38)
1/2
u
-
u
h
4.14e − 3
3.00e − 4
2.11e − 5
1.45e − 6
1.91 (1.25)
p
-
p
h
3.44e − 2
4.31e − 4
5.38e − 5
6.73e − 6
1.50 (1.25)
q
-
q
h
H
8.33e − 3
1.47e − 3
2.59e − 4
4.59e − 5
1.25 (1.25)
2/3
u
-
u
h
3.86e − 3
3.24e − 4
2.65e − 5
2.14e − 6
1.80 (1.00)
p
-
p
h
2.58e − 3
3.20e − 4
4.00e − 5
5.02e − 6
1.50 (1.17)
q
-
q
h
H
7.94e − 3
1.51e − 3
3.00e − 4
5.94e − 5
1.17 (1.17)Example 4.
Let source termf
(
x
)
=
x
-
1
/
4, which is singular at the origin and belongs to H
s
(
Ω
) for s
∈
[
0
,
1
/
4
). The analytic solution u
(
x
)
=
c
β
(
x
1
-
β
-
x
7
/
4
-
β
) only belongs to H
1
-
β
+
γ
(
Ω
) with γ
∈
[
β
/
2
,
1
/
2
) and c
β
=
Γ
(
3
/
4
)
/
Γ
(
11
/
4
-
β
). We split (1a) and (1b) into two equations (16a) and (16b) and (17a) and (17b). We are to solve (16a) and (16b) with right term f
r
(
x
)
=
x
-
1
/
4
-
2
c
β
/
Γ
(
1
+
β
)
x
β and analytic solution u
r
(
x
)
=
c
β
(
x
2
-
x
7
/
4
-
β
)
∈
H
7
/
4
-
β
+
ϵ
(
Ω
). The numerical results are shown in Tables 7 and 8.Table 7
Numerical result of Example4 on a uniform mesh of size h for k
=
1.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
7.65e − 2
1.91e − 2
4.79e − 3
1.20e − 3
1.00 (1.00)
p
-
p
h
1.61e − 2
3.15e − 3
5.76e − 4
1.0e2 − 4
1.22 (1.13)
q
-
q
h
H
2.53e − 2
5.78e − 3
1.24e − 3
2.63e − 4
1.10 (1.13)
1/2
u
-
u
h
7.38e − 2
1.84e − 2
4.60e − 3
1.15e − 3
1.00 (1.00)
p
-
p
h
1.38e − 2
2.70e − 3
4.91e − 4
8.76e − 5
1.22 (1.00)
q
-
q
h
H
2.97e − 2
8.02e − 3
2.04e − 3
5.13e − 4
0.98 (1.00)
2/3
u
-
u
h
7.27e − 2
1.81e − 2
4.53e − 3
1.13e − 3
1.00 (0.75)
p
-
p
h
1.26e − 2
2.49e − 3
4.51e − 4
8.03e − 5
1.22 (0.92)
q
-
q
h
H
3.02e − 2
9.51e − 3
2.71e − 3
7.63e − 4
0.88 (0.92)Table 8
Numerical result of Example4 on a uniform mesh of size h for k
=
2.
β
h
1/4
1/16
1/64
1/256
Order
1/4
u
-
u
h
5.50e − 3
3.99e − 4
2.80e − 5
1.92e − 6
1.91 (1.38)
p
-
p
h
1.14e − 2
2.05e − 3
3.62e − 4
6.40e − 5
1.25 (1.13)
q
-
q
h
H
1.81e − 2
3.86e − 3
8.12e − 4
1.71e − 4
1.12 (1.13)
1/2
u
-
u
h
3.95e − 3
3.63e − 4
3.27e − 5
2.92e − 6
1.73 (1.00)
p
-
p
h
3.91e − 2
6.88e − 4
1.22e − 4
2.15e − 5
1.25 (1.00)
q
-
q
h
H
1.14e − 2
2.84e − 3
7.09e − 4
1.77e − 4
1.00 (1.00)
2/3
u
-
u
h
2.36e − 3
2.65e − 4
2.96e − 5
3.31e − 6
1.58 (0.75)
p
-
p
h
2.48e − 3
4.33e − 4
7.65e − 5
1.38e − 5
1.25 (0.92)
q
-
q
h
H
7.12e − 3
1.99e − 3
5.60e − 4
1.58e − 4
0.92 (0.92)Remark 5.
From the above four examples, we can see that the error rates forq
-
q
h
H are identical to the predicted rates of convergence in Theorem 26. However, numerical convergence rates for p
-
p
h and u
-
u
h are much higher than the theoretical prediction of Theorem 26.



---

*Source: 9398265-2016-08-16.xml*