SUPERCONVERGENCE OF A DISCONTINUOUS ... - Semantic Scholar
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD FOR FRACTIONAL DIFFUSION AND WAVE EQUATIONS∗ KASSEM MUSTAPHA† AND WILLIAM MCLEAN‡ Abstract. We consider an initial-boundary value problem for ∂t u − ∂t−α ∇2 u = f (t), that is, for a fractional diffusion (−1 < α < 0) or wave (0 < α < 1) equation. A numerical solution is found by applying a piecewise-linear, discontinuous Galerkin (DG) method in time combined with a piecewiselinear, conforming finite element method in space. The time mesh is graded appropriately near t = 0, but the spatial mesh is quasiuniform. Previously, we proved that the error, measured in the spatial L2 -norm, is of order k2+α− + h2 `(k), uniformly in t, where k is the maximum time step, h is the maximum diameter of the spatial finite elements, α− = min(α, 0) ≤ 0 and `(k) = max(1, | log k|). Here, we prove convergence of order k3+2α− + h2 at each time level tn , for −1 < α < 1. Thus, if −1/2 < α < 1 then the DG solution is superconvergent, which generalizes a known result for the classical heat equation (i.e., the case α = 0). A simple postprocessing step employing Lagrange interpolation leads to superconvergence for any t. Numerical experiments indicate that our theoretical error bound is pessimistic if α < 0. Ignoring logarithmic factors, we observe that the error in the DG solution at t = tn , and after postprocessing at all t, is of order k3+α− + h2 for −1 < α < 1. Key words. finite elements, dual problem, postprocessing AMS subject classifications. 26A33, 35R09, 45K05, 47G20, 65M12, 65M15, 65M60
1. Introduction. In previous work [22, 30, 31, 32], we have studied discontinuous Galerkin (DG) methods for the time discretization of the abstract intial value problem u0 + Bα Au = f (t)
for 0 < t < T ,
with u(0) = u0 ,
(1.1)
where u0 = ∂u/∂t and Bα = ∂t−α ; more precisely, letting ωµ (t) = tµ−1 /Γ(µ) for µ > 0, the function Bα v is either a (Riemann–Liouville) fractional order derivative in time, ∂ Bα v(t) = ∂t
t
Z
ω1+α (t − s)v(s) ds if −1 < α < 0,
(1.2)
0
or a fractional order integral in time, Z Bα v(t) =
t
ωα (t − s)v(s) ds if 0 < α < 1. 0
In Section 2 we set out technical assumptions on the operator A, but for the present discussion we simply take Au = −∇2 u on a spatial domain Ω, and impose homogeneous Dirichlet boundary conditions on u. Problems of the form (1.1) arise in a variety of physical, biological and chemical applications [12, 18, 26, 27, 34, 38, 39, 40]. The case −1 < α < 0 describes slow or anomalous sub-diffusion and occurs, for example, in models of fractured or porous media, where the particle flux depends on the entire history of the density gradient ∇u. The case 0 < α < 1 describes wave propagation in viscoelastic materials [10, 17, 35]. ∗ We
thank KFUPM for supporting this research as part of the project SB101020. Department of Mathematics and Statistics, KFUPM, Dhahran 31261, Saudi Arabia ([email protected]). ‡ School of Mathematics and Statistics, The University of New South Wales, Sydney 2052, Australia ([email protected]). †
1
2
KASSEM MUSTAPHA AND WILLIAM MCLEAN
In the limit as α → 0, the evolution equation in (1.1) becomes u0 + Au = f , which is just the classical heat equation, and Eriksson et al. [9] studied the convergence of the DG solution U (t) ≈ u(t) in this case. For a maximum time step k, and using discontinuous piecewise polynomials of degree at most q − 1 in t, with no spatial discretization, they proved for 0 ≤ t ≤ T an optimal convergence rate Z t kf (q) (s)k ds , kU (t) − u(t)k ≤ Ck q ku0 kq + ku(q) (0)k + kf (q−1) (0)k + 0
where kvk is the norm of v in L2 (Ω) and where kvkq = kAq/2 vk for v in D(Aq/2 ), the domain of the unbounded operator Aq/2 . In addition, they proved that the DG solution is superconvergent at the nth time level tn , satisfying an error bound Z tn (q) − 2q−1 (q) kf (s)kq−1 ds , kU (tn ) − u(tn )k ≤ Ck ku0 k2q−1 + ku (0)kq−1 + 0
U (t) denotes the limit from the left. Ericksson et al. were also where U (t− n ) = limt→t− n able to prove that a convergence rate faster then O(k q ) holds under less restrictive spatial regularity requirements on the solution u. Our aim is to establish superconvergence results for the fractional-order problem (1.1), restricting our attention to the piecewise-linear DG method (q = 2). We believe that this time-stepping scheme is the first to achieve better than second-order accuracy. As well as nodal superconvergence of the DG solution we show that a postprocessed solution is superconvergent uniformly in t. Many authors have studied numerical methods for (1.1). In the case 0 < α < 1, Sanz-Serna [36] proposed a convolution quadrature scheme, and subsequently Cuesta, Lubich and Palencia [4, 5, 3] developed this approach to obtain an O(k 2 ) method as well as a fast implementation [37]. McLean and Thom´ee [23] combined finite differences and quadrature in time, with finite elements in space. In the case −1 < α < 0, Langlands and Henry [13] introduced an implicit Euler scheme involving the Gr¨ unwald–Letnikov fractional derivative and spatial finite differences with step size h, and observed O(k 1/2 + h2 ) convergence in the case α = −1/2. Yuste and Acedo [43] treated an explicit Euler scheme and showed O(k + h2 ) convergence. Zhuang, Liu, Anh, Turner et al. [2, 14, 44, 45] developed another class of O(k + h2 ) finite difference methods, and Yuste [42] presented an O(k 2 + h2 ) method. Cui [6] and Chen et al. [1] studied O(k + h4 ) schemes, and Cui [7, 8] analysed an O(k min(1−α,2+α) + h4 ) ADI scheme on a rectangular spatial domain; see also Wang and Wang [41] and Zhang and Sun [33]. For another type of finite difference scheme [28, 29], the error is O(k 2+α + h2 ), and recently Jin et al. [11] proved optimal error bounds for two semidiscrete finite element methods. Some of these works employ an alternative formulation of (1.1) using the Caputo fractional derivative. In practice, the higher order derivatives of u are typically singular [19, 21] as t → 0, so formally high order methods [1, 2, 6, 7, 8, 14, 33, 41, 43, 44, 45] can fail to achieve fast convergence. We have analysed several methods that allow for the singular behaviour of u by employing non-uniform time steps [21, 25, 28, 29, 32]. Another approach, that yields a parallel in time algorithm with spectral accuracy even for problems with low regularity, is to approximate u via the Laplace inversion formula [15, 16, 24]. To minimise the need for handling separately the cases α < 0 and α > 0, it is convenient to write α+ = max(α, 0) ≥ 0 and α− = min(α, 0) ≤ 0 for the positive and
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
3
negative parts of α, respectively. To prove our superconvergent nodal error bounds, we assume that there exist positive constants M and σ such that kAu0 k + kAu(t)k ≤ M
and kAu0 (t)k + tkAu00 (t)k ≤ M tσ−1 ,
(1.3)
as well as tkA2 u0 (t)k + t2 kA2 u00 (t)k ≤ M tσ−α− −1 ,
(1.4)
for 0 < t ≤ T . For instance [19, 21], if f ≡ 0 and u0 ∈ D(A2 ), then (1.3) and (1.4) hold with M = CkAu0 k and σ = 1 + α− . Section 2 sets out our notation and assumptions, and recalls some tools and results from earlier work [31]. In Section 3, we introduce the homogeneous dual problem, −z 0 + Bα∗ Az = 0
for 0 < t < T ,
with z(T ) = zT ,
(1.5)
for a given terminal value zT , and represent the nodal error U (t− n ) − u(tn ) in terms of z(t) and its DG approximation Z(t). We allow a class of non-uniform meshes, specified in Section 4, where we prove in Theorem 4.3 that the nodal error is O(k 3+2α− ). Our method of analysis allows us to handle the two cases −1 < α < 0 and 0 < α < 1 together, but the former presents additional technical difficulties in some places. In an earlier paper [30, Theorem 4.1], we estimated the nodal error for the case 0 < α < 1 in a different way that yields a bound of order k 2+α . (Although we claimed O(k 3 ) convergence, the first line of [30, Corollary 4.2] contains an error.) In Section 5 we construct, via a simple interpolation scheme, a postprocessed solution U ] whose error is O(k 3+2α− ) for all t, not just at the nodal values, assuming that ku0 (t)k + t2 ku000 (t)k ≤ M tσ
]
−1
for 0 < t ≤ T .
(1.6)
For instance [19], if u0 ∈ D(A) and f ≡ 0, then (1.6) holds with σ ] = 1 + α and M = kAu0 k. Section 6 introduces a fully discrete scheme by applying a continuous piecewise-linear, finite element method for the spatial discretization. Thus, the fully discrete solution is continuous in space but discontinuous in time. We show that the error bound is as for the semidiscrete method but with an extra term of order h2 . Finally, we present some numerical examples in Section 7, which indicate that our error bounds are pessimistic, at least in some cases. We observe that the nodal error from the time discretization is O(k 3+α− ), which is better than our theoretical estimate by a factor k α− . The same is true for the postprocessed solution, uniformly in t. 2. Preliminaries. 2.1. Assumptions on the spatial operator. We assume as in earlier work [9, 31] that the self-adjoint linear operator A has a complete eigensystem in a real Hilbert space H, say Aφj = λj φj for j = 1, 2, 3, . . . , and that A is strictly positive-definite with the eigenvalues ordered so that 0 < λ1 ≤ λ2 ≤ · · · . (Strict positive definiteness is not essential, but allowing λ1 = 0 would result in some technical complications [30, Section 5] that we prefer to avoid here.) We denote p the inner product of u and v in H by hu, vi and the corresponding norm by kuk = hu, ui. Associated with the linear operator A is a bilinear form, denoted by the same symbol: A(u, v) =
∞ X m=1
λm hu, φm ihφm , vi for u, v ∈ D(A1/2 ).
4
KASSEM MUSTAPHA AND WILLIAM MCLEAN
These assumptions hold, in particular, if A = −∇2 subject to homogenous Dirichlet boundary conditions on R a bounded domain Ω, because A has a compact inverse on H = L2 (Ω) and A(u, v) = Ω ∇u · ∇v dx. 2.2. The discontinuous Galerkin time discretization. Fixing a time interval [0, T ], we introduce a mesh for the time discretization, 0 = t0 < t1 < · · · < tN = T,
(2.1)
with kn = tn − tn−1 and In = (tn−1 , tn ) for 1 ≤ n ≤ N , and a maximum time step k = max1≤n≤N kn . Let Pr denote the space Sn of polynomials of degree at most r with coefficients in D(A1/2 ), and let Jn = j=1 Ij = [0, tn ] \ {t0 , t1 , . . . , tn }, with J = JN . Our trial space W consists of the piecewise-linear functions U : J → D(A1/2 ) with U |In ∈ P1 for 1 ≤ n ≤ N . We treat U as undefined at each time level tn , and write n U− = U (t− n ),
n U+ = U (t+ n ),
n n [U ]n = U+ − U− .
(2.2)
For r ∈ {0, 1, 2, . . .} we let C r (J, H) denote the space of functions v : J → H such that the restriction v|In extends to an r-times continuously differentiable function on the closed interval [tn−1 , tn ], for 1 ≤ n ≤ N . In other words, v is a piecewise C r function with respect to the time levels tn . If v ∈ C 1 (J, H), then its fractional derivative (1.2) admits the representation [31] Bα v(t) =
0 ω1+α (t)v+
+
n−1 X
j
ω1+α (t − tj )[v] +
j=1
n Z X j=1
min(tj ,t)
ω1+α (t − s)v 0 (s) ds (2.3)
tj−1
for t ∈ In and −1 < α < 0. Thus, Bα v(t) is left-continuous at t = tn−1 but has a n−1 weak singularity (t − tn−1 )α as t → t+ 6= 0. However, if 0 < α < 1 then n−1 if [v] Bα v(t) is continuous for 0 ≤ t ≤ T . For −1 < α < 1, the piecewise-linear DG time stepping procedure determines 0 U ∈ W by setting U− = u0 and requiring [30, 31] n−1 n−1 i , X+ hU+
Z +
0 hU (t), X(t)i + A Bα U (t), X(t) dt
In
=
n−1 n−1 i , X+ hU−
Z hf (t), X(t)i dt, (2.4)
+ In
for 1 ≤ n ≤ N and for every test function X ∈ P1 . The nonlocal nature of the operator Bα means that at each time step we must compute a sum involving all previous times levels, but this sum can be evaluated via a fast algorithm [20]. 2.3. Galerkin orthogonality and stability. For v ∈ C 1 J, D(A1/2 ) and w ∈ C J, D(A1/2 ) , we define the global bilinear form GN (v, w) =
0 0 hv+ , w+ i
+
N −1 X
h[v]
n
n , w+ i
n=1
N Z X 0 + hv , wi + A(Bα v, w) dt. n=1
(2.5)
In
Summing the equations (2.4) gives 0 0 GN (U, X) = hU− , X+ i+
Z
tN
hf (t), X(t)i dt for all X ∈ W, 0
(2.6)
5
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
and conversely, (2.6) implies that U satisfies (2.4) for 1 ≤ n ≤ N . Since [u]n = 0, Z tN 0 GN (u, X) = hu0 , X+ i + hf (t), X(t)i dt, (2.7) 0 0 and thus, assuming U− = u0 , the error has the Galerkin orthogonality property
GN (U − u, X) = 0
for all X ∈ W.
(2.8)
The DG method is unconditionally stable. Indeed, with the notation kU kI = sup kU (t)k
for any I ⊆ [0, T ],
t∈I
the following estimate holds. 0 Theorem 2.1. Given U− ∈ H and f ∈ L1 (0, T ); H , there exists a unique U ∈ W satisfying (2.4) for n = 1, 2, . . . , N . Furthermore, U (t) ∈ D(A) (t > 0), and Z tn 0 0 2 kU kJn ≤ 8 hU− , U+ i + hf (t), U (t)i dt for 1 ≤ n ≤ N . 0
Proof. Since hV 0 , V i = 12 (d/dt)kV k2 and 2GN (V, V ) ≥ kV+0 k2 + kV−N k2 +
RT 0
N −1 X
A(Bα V, V ) dt ≥ 0 we find that k[V ]n k2
for all V ∈ W,
(2.9)
n=1
implying the stated estimate; for details, see [30, Theorem 2.1], [31, Theorem 1]. 2.4. A discontinuous quasi-interpolant. The conditions Z − Π− v(t− ) = v(t ) and v(t) − Π− v(t)] dt = 0 n n
(2.10)
In
determine a unique projection operator Π− : C(J, H) → W; explicitly, Π− v(t) := v(t− n)+
¯n v(t− n)−v (t − tn ) kn /2
for t ∈ In ,
R where v¯n = kn−1 In v(t) dt denotes the mean value of v over In . The interpolation error admits the integral representations [30, Equation (3.8)] Z tn Z tn − t Π− v(t) − v(t) = v 0 (s) ds − 2 (s − tn−1 )v 0 (s) ds 2 k t In n (2.11) Z tn Z t − t n 00 2 00 = (t − s)v (s) ds + (s − tn−1 ) v (s) ds, for t ∈ In . kn2 t In Likewise, the conditions + Π+ v(t+ n−1 ) = v(tn−1 )
Z and
v(t) − Π+ v(t)] dt = 0,
In
determine a unique projector Π+ : C(J, H) → W, with Π+ v(t) := v(t+ n−1 ) +
v¯n − v(t+ n−1 ) (t − tn−1 ) kn /2
for t ∈ In ,
(2.12)
6
KASSEM MUSTAPHA AND WILLIAM MCLEAN
and t
Z t − tn−1 v (s) ds + 2 Π v(t) − v(t) = − (tn − s)v 0 (s) ds kn2 tn−1 In Z Z t t − tn−1 = (s − tn−1 )v 00 (s) ds + (tn − s)2 v 00 (s) ds. 2 k tn−1 In n Z
+
0
(2.13)
Thus, for t ∈ In (Π+ v − v)0 (t) = −v 0 (t) +
2 kn2
Z
(tn − s)v 0 (s) ds
(2.14)
In
and short calculations lead to the error bounds Z kv (r) (t)k dt ≤ (4−r)knr kv (r) kIn kΠ± v −vkIn ≤ (4−r)knr−1
for r ∈ {1, 2}, (2.15)
In
and the stability estimates Z
± 0
(Π ) v ≤ 2 In kn
kΠ± vkIn ≤ 3kvkIn ,
kv 0 (t)k dt,
± n
[Π v] ≤
Z
kv 0 (t)k dt.
In
In
(2.16) 3. Dual problem. 3.1. Properties of the adjoint operator. The adjoint operator appearing in the dual problem (1.5) should satisfy, for appropriate u and v, the identity Z
T
Z
T
hv, Bα wi dt = 0
hBα∗ v, wi dt,
(3.1)
0
and the next lemma establishes an explicit representation of Bα∗ . Lemma 3.1. The identity (3.1) holds in the following cases. 1. If −1 < α < 0 and v, w ∈ C 1 (J, H), with Bα∗ w(t) = −
∂ ∂t
T
Z
ω1+α (s − t)w(s) ds
for t ∈ J.
t
2. If 0 < α < 1 and v, w ∈ C(J, H), with Bα∗ w(t)
Z
T
ωα (s − t)w(s) ds
=
for 0 ≤ t ≤ T .
t
Proof. In case 1, we see from the representation (2.3) that Z
T
hBα v, wi dt = 0
where, letting B n,j =
N Z X n=1
R In
hBα v, wi dt = S1 + S2 + S3 ,
In
ω1+α (t − tj )w(t) dt and
Z
Z
min(tj ,t)
Dnj = In
tj−1
ω1+α (t − s)hv 0 (s), w(t)i ds dt,
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
7
we define S1 =
N X
0 n,0 v+ , B ,
S2 =
n=1
N n−1 X X
[v]j , B n,j ,
S3 =
n=2 j=1
N X n X
Dnj .
n=1 j=1
By reversing the order of integration, integrating by parts and then interchanging variables, we find that for 1 ≤ j ≤ n − 1, Z Z ∂ j j−1 ω1+α (s − t)w(s) ds dt, Dnj = hv− , B n,j i − hv+ , B n,j−1 i − v(t), ∂t In Ij whereas n−1 Dnn = −hv+ , B n,n−1 i −
Z
∂ ∂t
v(t),
In
tn
Z
ω1+α (s − t)w(s) ds dt.
t
Thus, after interchanging the order of summation for the double integrals, N Z N n−1 N X j X X X j−1 n−1 n,j n,j−1 n,n−1 hv, Bα∗ wi dt, hv− , B i − hv+ , B i − i+ S3 = − hv+ , B n=2 j=1
n=1
j=1
Ij
that is, S3 =
N n−1 X X
j 0 hv− , B n,j i − hv+ , B 1,0 i −
n=2 j=1
Z N n−1 X X j hv+ , B n,j i −
hv, Bα∗ wi dt
0
n=2 j=0
Z
T
T
hv, Bα∗ wi dt,
= −S1 − S2 − 0
so (3.1) holds. In the case 0 < α < 1, we simply reverse the order of integration. The adjoint operator admits a representation analogous to (2.3). Lemma 3.2. If −1 < α < 0, then Bα∗ w(t) equals ω1+α (tn −
N t)w−
−
N −1 X
j
ω1+α (tj − t)[w] −
N Z X j=n
j=n
tj
ω1+α (s − t)w0 (s) ds
max(tj−1 ,t)
Bα∗ w(t) is n t− n if [w]
1
for w ∈ C (J, H) and t ∈ In . Thus, right-continuous at t = tn but possesses a weak singularity (tn − t)α as t → 6= 0. Proof. If t ∈ In and j ≥ n + 1, then, integrating by parts, Z Z j j−1 ω1+α (s−t)w(s) ds = ω2+α (tj −t)w− −ω2+α (tj−1 −t)w+ − ω2+α (s−t)w0 (s) ds Ij
Ij
and Z
tn
ω1+α (s − t)w(s) ds = ω2+α (tn −
n t)w−
Z
tn
−
t
ω2+α (s − t)w0 (s) ds.
t
Differentiating these expressions with respect to t, we see from part 1 of Lemma 3.1 that Bα∗ w(t) equals N X j=n
j ω1+α (tj − t)w−
−
N X j=n+1
j−1 ω1+α (tj−1 − t)w+
−
N Z X j=n
tj
max(tj−1 ,t)
and the result follows after shifting the index in the second sum.
ω1+α (s − t)w0 (s) ds,
8
KASSEM MUSTAPHA AND WILLIAM MCLEAN
3.2. Representation of the nodal error. Integration by parts in (2.5), together with the identity (3.1), shows that for all v, w ∈ C 1 (J, H), GN (v, w) =
N N hv− , w− i
−
N −1 X
n hv− , [w]n i
N Z X + −hv, w0 i + A(v, Bα∗ w) dt.
n=1
n=1
(3.2)
In
Since −hv, z 0 i + A(v, Bα∗ z) = hv, −z 0 + Bα∗ Azi = 0, the solution z of the dual problem (1.5) satisfies N GN (v, z) = hv− , zT i for all v ∈ C J, D(A1/2 ) . (3.3) We therefore define the DG solution Z ∈ W of (1.5) by N GN (V, Z) = hV−N , Z+ i for all V ∈ W,
(3.4)
N with Z+ = zT , and deduce the Galerkin orthogonality property
GN (V, Z − z) = 0
for all V ∈ W.
(3.5)
The following representation is the basis for our analysis of the nodal error. Theorem 3.3. If u and z are the solutions of the initial-value problem (1.1) and of the dual problem (1.5), and if U and Z are the corresponding DG solutions, then N hU− − u(tN ), zT i = GN (u − Π− u, Z − z)
for every zT ∈ H.
Proof. Taking V = U in (3.4) and v = u in (3.3) gives N N hU− − u(tN ), zT i = hU− , zT i − hu(tN ), zT i = GN (U, Z) − GN (u, z) = GN (u, Z − z),
where the last step used the Galerkin orthogonality property (2.8) of U , with X = Z. Now use the Galerkin orthogonality property (3.5) of Z, with V = Π− u. 3.3. Error in the DG solution of the dual problem. We will use the following regularity estimates. Lemma 3.4. For −1 < α < 1 and 0 < t < T , the solution z of the dual problem (1.5) satisfies kA−1 z 0 (t)k + (T − t)kA−1 z 00 (t)k ≤ C(T − t)α kzT k and (T − t)1+α kAz(t)k + kz(t)k + (T − t)kz 0 (t)k ≤ CkzT k. Proof. Define the time reversal operator Rv(t) = v(T − t). Since R∂t = −∂t R and RBα∗ = Bα R, we deduce from (1.5) that the function v = RA−1 z satisfies v 0 + Bα Av = 0
for 0 < t < T ,
with v(0) = A−1 zT .
Known results for −1 < α < 0 [19, Theorem 4.2] and 0 < α < 1 [21, Theorem 2.1 with r = 2 and µ = −2] give kv 0 (t)k + tkv 00 (t)k ≤ Ctα kAv(0)k = Ctα kzT k, implying the first estimate. Similarly [19, Theorem 4.1], the function w = Rz satisfies t1+α kAw(t)k + kw(t)k + tkw0 (t)k ≤ Ckw(0)k = CkzT k,
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
9
implying the second estimate. To investigate the DG error for the dual problem, we make the splitting A−1 (Z − z) = ζ + Θ
where ζ = A−1 (Π+ z − z) and Θ = A−1 (Z − Π+ z) ∈ W. (3.6)
Lemma 3.5. The function ζ in (3.6) satisfies α
kζkJ ≤ CtN+ k 1+α− kzT k . Proof. By (2.15) and Lemma 3.4, kζkJ is bounded by Z Z k(I − Π+ )A−1 zkJ ≤ 3 max kA−1 z 0 (t)k dt ≤ CkzT k max (tN − t)α dt. 1≤n≤N
1≤n≤N
In
In
If −1 < α < 0, then Z (1 + α) (tN − t)α dt = (tN − tn−1 )1+α − (tN − tn )1+α ≤ kn1+α , In
R whereas if 0 < α < 1, then In (tN − t)α dt ≤ kn tα N. Lemma 3.6. The function Θ in (3.6) satisfies kΘk2J
Z ≤ 8
T
0
hΘ(t), Bα∗ Aζi dt .
Proof. By (3.5), GN (V, ζ + Θ) = G(A−1 V, Z − z) = 0 for all V ∈ W, where we used the identity GN (v, A−1 w) = GN (A−1 v, w) and the fact that A−1 V ∈ W. Thus, GN (V, Θ) = −G(V, ζ)
for all V ∈ W.
n The first property in (2.12) means that ζ+ = 0 for 0 ≤ n ≤ N − 1, so the formula (3.2) shows
GN (V, ζ) =
N X
n hV−n , ζ− i+
n=1
N Z X −hV, ζ 0 i + A(V, Bα∗ ζ) dt, n=1
In
R R n n and integration by parts gives In hV, ζ 0 i dt = hV−n , ζ− i − In hV 0 , ζi dt = hV−n , ζ− i, 0 where, in the last step, we used the second property in (2.12) and the fact that V is constant on In . Thus, if we define g = −ABα∗ ζ = −Bα∗ Aζ then Z GN (V, Θ) = − 0
T
A(V, Bα∗ ζ) dt
Z
T
hV, gi dt
=
for all V ∈ W,
0
which means that Θ ∈ W is the DG solution of −θ0 + Bα∗ Aθ = g(t) for 0 < t < T , with θ(T ) = 0. The desired estimate follows by the stability of Θ, which we can prove by applying Theorem 2.1 to RΘ(t) = Θ(T − t). Lemma 3.7. If −1 < α < 1 then Z T α ∗ hV, Bα Aζi dt ≤ CtN+ k 1+α− `(tN /kN )kV kJ kzT k for all V ∈ W. 0
10
KASSEM MUSTAPHA AND WILLIAM MCLEAN
n−1 Proof. Suppose first that −1 < α < 0. Since Bα V = (B1+α V )0 and ζ+ = 0, we see using (3.1) and integrating by parts that N Z X
T
Z
hV, Bα∗ Aζi dt =
0
h(B1+α V )0 , Aζi dt =
In
n=1
N Z X
h∆n , Aζ 0 i dt,
In
n=1
where, for t ∈ In , ∆n (t) = B1+α V (tn ) − B1+α V (t) Z Z t [ω1+α (tn − s) − ω1+α (t − s)]V (s) ds + =
tn
ω1+α (tn − s)V (s) ds.
t
0
The function ω1+α is monotone decreasing whereas ω2+α is monotone increasing, so t
Z tn ω1+α (tn − s) ds k∆ (t)k ≤ kV kJ ω1+α (t − s) − ω1+α (tn − s) ds + 0 t = kV kJ ω2+α (t) − ω2+α (tn ) + 2ω2+α (tn − t) ≤ 2kV kJ ω2+α (tn − t) Z
n
R T and the Cauchy–Schwarz inequality shows that 0 hV, Bα∗ Aζi dt is bounded by 2kV kJ
NX −1
Z
ω2+α (tN − t)kAζ k dt .
Z
0
0
kAζ k dt +
ω2+α (kn ) In
n=1
IN
Since Aζ = Π+ z − z, the representation (2.14) and Lemma 3.4 imply N −1 X
Z
0
kAζ k dt ≤ C
ω2+α (kn ) In
n=1
=C
N −1 X
kn1+α
Z
n=1
In
N −1 X
Z
kn1+α
2 kz (t)k + 2 kn 0
1+
≤ Ck 1+α kzT k
Z
(tn − s)kz (s)k ds dt 0
In
In
n=1
Z
tN −1
2(tn − t) kz 0 (t)k dt kn
(T − t)−1 dt = Ck 1+α log(tN /kN )
0
and Z IN
Z 2 0 0 (T −s)kz (s)k ds dt ω2+α (tN −t)kAζ k dt ≤ C (T −t) kz (t)k+ 2 kN IN IN Z α 2kN 1+α =C (T − t)α + (T − t)kz 0 (t)k dt ≤ CkzT kkN . 2+α IN Z
0
1+α
The desired estimate follows at once. Now let 0 < α < 1. By part 2 of Lemma 3.1, Z
0
T
Z hV, Bα∗ Aζi dt ≤
T
Z kV (t)k
0
Z ≤ kV kJ
T
ωα (s − t)kAζ(s)k ds dt t
Z kAζ(s)k
0
T
s α
Z
ωα (s − t) dt ds ≤ CT kV kJ 0
T
kAζ(t)k dt. 0
11
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
The estimates (2.15) and (2.16) imply that T
Z
kAζ(t)k dt ≤ 0
N X
kn kAζkIn ≤ 4kN kzkIN + 3
n=1
N −1 X
Z kn
n=1
kz 0 (t)k dt,
In
and we know from Lemma 3.4 that kzkIN ≤ CkzT k and N −1 X
Z
Z
0
tN −1
kz (t)k dt ≤ CkkzT k
kn In
n=1
(tN − t)−1 dt = CkkzT klog(tN /kN ).
0
Hence, we arrive at the following error estimate for the dual problem. Theorem 3.8. Let z denote the solution of the dual problem (1.5), and let Z denote the DG solution defined by (3.4). Then, for −1 < α < 1, α
kA−1 (Z − z)kJ ≤ CtN+ k 1+α− `(tN /kN )kzT k. Proof. The splitting (3.6) implies that kA−1 (Z − z)kJ ≤ kζkJ + kΘkJ , and we estimate these two terms using Lemmas 3.5, 3.6 and 3.7. 4. Nodal superconvergence. With the help of Theorems 3.3 and 3.8, we are n ≈ u(tn ). Define now able to estimate the error in the approximation U− (u) = D1 + E1 + max kj Dj + 2≤j≤n
n X
2+α−
kj
Ej ,
(4.1)
j=2
where Z D1 =
kAu0 (t)k dt
Z and E1 =
t1+α− kA2 u0 (t)k dt,
(4.2)
I1
I1
with Z Dn =
kAu00 (t)k dt
Z and En =
In
kA2 u00 (t)k dt
for 2 ≤ n ≤ N .
(4.3)
In
For the rest of the paper, we assume that k1 ≤ k2 ≤ · · · ≤ kN .
(4.4)
Theorem 4.1. Let u be the solution of the initial value problem (1.1) and let U be the DG solution satisfying (2.4). Then, for 1 ≤ n ≤ N , n α+ 1+α− + kU− − u(tn )k ≤ Ctα `(tn /kn ) (u). n (1 + tn )k
(4.5)
Proof. Put η = u − Π− u and define Z Z n 0 n δ1 = hη, (z − Z) i dt and δ2 = hBα Aη, z − Zi dt. In
Since
n η−
In
= 0 for 1 ≤ n ≤ N , we see from Theorem 3.3, (3.2) and Lemma 3.1 that N hU− − u(tN ), zT i = GN (η, Z − z) =
N X n=1
δ1n + δ2n .
12
KASSEM MUSTAPHA AND WILLIAM MCLEAN
Since Z 0 is constant on In , the second property of Π− in (2.10) gives Z Z Z Z t
n 0 0 0 Aη(t), A−1 z 00 (s) ds dt, δ1 = hη, z i dt = η(t), z (t) − z (tn−1 ) dt = In
In
tn−1
In
and therefore, using Lemma 3.4, N −1 X
|δ1n |
N −1 X
≤ kAηkJ
n=1
Z
whereas N −1 X
Z kn
R =
z (t)k dt ≤ CkAηkJN kzT k
N −1 X
In
n=1
|δ1N |
kA
kn
−1 00
Z
(tN −t)α−1 dt,
kn
n=1
In
R z i dt ≤ CkAηkIN kzT k IN (tN − t)α dt. Here,
−1 0
IN
hAη, A
(tN − t)α−1 dt
In
n=1
−α−
≤ k 1+α− kN
−α−
tN −1
Z
(tN − t)α−1 dt =
0
k 1+α− kN α
α
+ 1+α− α (tα , N − kN ) ≤ CtN k
R α 1+α and likewise IN (tN − t)α dt = kN /(1 + α) ≤ CtN+ k 1+α− . By (2.15), kAηkI1 ≤ 3D1 and kAηkIn ≤ 2kn Dn for 2 ≤ n ≤ N , so N X
α |δ1n | ≤ CtN+ k 1+α− kzT k D1 + max kn Dn for −1 < α < 1. 2≤n≤N
n=1
(4.6)
Turning to δ2n , if −1 < α < 0, then [31, Lemma 2] N Z X n δ2 = n=1
tN 2
−1
hBα A η, A
0
N X −1 2+α (Z − z)i dt ≤ CkA (Z − z)kJ E1 + kn En , n=2
(4.7) but if 0 < α < 1 then |δ2n | is bounded by Z Z −1 2 −1 kA (Z − z)kIn kBα A η(t)k dt ≤ kA (Z − z)kIn In
In
Z
t
ωα (t − s)kA2 η(s)k ds,
0
so, after summing over n and reversing the order of integration, N X
−1 |δ2n | ≤ Ctα (Z − z)kJ N kA
Z
tN
kA2 η(t)k dt.
0
n=1
The integral representation (2.11) implies that Z Z Z Z t1 t1 − t skA2 u0 (s)k ds dt kA2 η(t)k ≤ kA2 u0 (s)k ds + 2 2 k1 I1 I1 I t Z s Z1 Z s 2 0 = kA u (s)k dt + 2 2(t1 − t) dt ds = 2E1 , k1 I1 I1 0 Rt PN PN and by (2.15), t1N kA2 η(t)k dt ≤ n=2 kn kA2 ηkIn ≤ n=2 2kn2 Ej . Combine this with (4.7) and applying Theorem 3.8, N X n=1
N X 2α |δ2n | ≤ CtN + k 1+α− `(tN /kN )kzT k E1 + kn2+α− En for −1 < α < 1. (4.8) n=2
13
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
Since zT ∈ H is arbitrary, the desired estimate follows from (4.6) and (4.8). To estimate the convergence rate at the nodes, we introduce some assumptions about the behaviour of the time steps, namely that, for some fixed γ ≥ 1, kn ≤ Cγ k min(1, t1−1/γ ) n
and tn ≤ Cγ tn−1
for 2 ≤ n ≤ N ,
(4.9)
with cγ k γ ≤ k1 ≤ Cγ k γ .
(4.10)
For example, these assumptions are satisfied if we put tn = (n/N )γ T
for 0 ≤ n ≤ N .
(4.11)
Recall that `(t) = max(1, | log t|). Lemma 4.2. Assume that u satisfies (1.3) and (1.4), and that the time mesh satisfies (4.9) and (4.10). Then, with γ ∗ = (2 + α− )/σ and for 1 ≤ n ≤ N , γσ 1 ≤ γ < γ∗, k , (u) ≤ CT M × k 2+α− `(tn /k1 ), γ = γ ∗ , 2+α− k , γ > γ∗. Rk Proof. The stated assumptions imply that D1 + E1 ≤ CM 0 1 tσ−1 dt ≤ CM k1σ ≤ CM k γσ , and, for 2 ≤ j ≤ n, ( Z k γσ , 1 ≤ γ < 2/σ, σ−2 2 σ−2 kj Dj ≤ CM kj t dt ≤ CM kj tj ≤ CM × σ−2/γ 2 k , γ ≥ 2/σ, tn Ij γσ 1 ≤ γ ≤ γ∗, k , 2+α − ≤ CM × k , γ ∗ ≤ γ ≤ 2/σ, σ−(2+α− )/γ 2+α− , γ ≥ 2/σ. tn k Similarly, n X j=2
2+α kj − Ej
≤M
n X
2+α kj −
j=2
≤ CM k 2+α−
Z t
σ−3−α−
dt ≤ CM k
2+α−
Z
Ij
tn
tσ−1−(2+α− )/γ dt
k1
γσ−(2+α− ) , 1 < γ < γ∗, k × log(tn /k1 ), γ = γ∗, σ−(2+α− )/γ tn , γ > γ∗.
We can now state our main result on nodal superconvergence. Theorem 4.3. Assume that the solution u of the initial value problem (1.1) satisfies (1.3) and (1.4), and that the time mesh satisfies (4.9) and (4.10). Then, with γ ∗ = (2 + α− )/σ and for 1 ≤ n ≤ N , the DG method (2.4) satisfies the nodal error bound 1+α− +γσ `(tn /kn ), 1 ≤ γ < γ ∗ , k n 3+2α − kU− − u(tn )k ≤ CT M × k `(tn /kn )2 , γ = γ∗, 3+2α− k `(tn /kn ), γ > γ∗ . Proof. The error bound follows at once from Theorem 4.1 and Lemma 4.2.
14
KASSEM MUSTAPHA AND WILLIAM MCLEAN
5. Postprocessing. We can postprocess the DG solution U to obtain a globally superconvergent solution U ] using simple Lagrange interpolation, as follows. Given a piecewise continuous function v : J → H, define Lv : J → H by linear interpolation on the first two subintervals, n−1 n Lv(t) = kn−1 [(tn − t)v− + (t − tn−1 )v− ]
for t ∈ In and n ∈ {1, 2},
(5.1)
and backward quadratic interpolation on the remaining subintervals, (t − tn−1 )(t − tn ) n−2 (t − tn−2 )(t − tn ) n−1 (t − tn−2 )(t − tn−1 ) n v − v− + v− kn−1 (kn−1 + kn ) − kn−1 kn (kn−1 + kn )kn (5.2) n for t ∈ In and n ≥ 3. Thus, (Lv)(tn ) = v− for 0 ≤ n ≤ N , and we define the postprocessed solution by Lv(t) =
U ] = LU.
(5.3)
The interpolant of the exact solution satisfies the following error bound. Lemma 5.1. Assume that there exist positive constants M and σ ] such that the regularity estimate (1.6) holds, and put γ ] = 3/σ ] . If the time mesh satisfies (4.9)– (4.10), then ( ] k γσ , 1 ≤ γ < γ], ku − LukJ ≤ CM × σ ] −3/γ 3 T k , γ ≥ γ]. Proof. If n ∈ {1, 2} and t ∈ In , then Z
t
u0 (s) ds −
(u − Lu)(t) = tn−1
t kn
Z
tn
u0 (s) ds,
tn−1
and thus, noting that t2 ≤ Cγ t1 , ku − LukIn ≤
Z tn ] ] tn 1+ ku0 (s)k ds ≤ CM tσ1 ≤ CM k γσ . kn tn−1
If n ≥ 3 and t ∈ In , then we can write the interpolation error in terms of a divided difference, (u − Lu)(t) = u[tn−2 , tn−1 , t, tn ](t − tn−2 )(t − tn−1 )(t − tn ), so ]
1 ku − LukIn ≤ 41 kn2 (kn−1 + kn ) 3! ku000 k[tn−2 ,tn ] ≤ CM kn2 (kn−1 + kn )tσn
−3
,
where, in the final step, we used (4.9). If 1 ≤ γ < 3/σ ] then, again using (4.9), ]
kn2 (kn−1 + kn )tσn
−3
]
]
]
≤ C(kt1−1/γ )γσ kn3−γσ tσn n
−3
]
]
]
= Ck γσ (kn /tn )3−γσ ≤ Ck γσ ,
but for γ ≥ 3/σ, ]
kn2 (kn−1 + kn )tσn
−3
]
≤ C(kt1−1/γ )3 tσn n
−3
]
≤ Ck 3 tσn
−3/γ
≤ CT σ
]
−3/γ 3
k .
Now consider the stability of the interpolation operator L. We see from (5.1) that 0 1 1 2 kLvkI1 ≤ max |v− |, |v− | and kLvkI2 ≤ max |v− |, |v− | .
15
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
A similar estimate holds for the subsequent subintervals provided the mesh satisfies the local quasi-uniformity condition kn ≤ Λkn−1
for 3 ≤ n ≤ N .
(5.4)
For example, our standard mesh (4.11) satisfies this condition with Λ = 2γ − 1. Lemma 5.2. If (5.4) holds, then j kLvkIn ≤ 2 + 45 Λ max |v− | for 3 ≤ n ≤ N . n−2≤j≤n
Proof. The estimate follows from (5.2) because, for t ∈ In and n ≥ 2, 1 2 1 kn |(t − tn−1 )(t − tn )| 4 kn ≤ ≤ 4 ≤ 14 Λ, kn−1 (kn−1 + kn ) kn−1 (kn−1 + kn ) kn−1 |(t − tn−2 )(t − tn )| (kn−1 + kn )kn kn ≤ =1+ ≤ 1 + Λ, kn−1 kn kn−1 kn kn−1 (kn−1 + kn )kn |(t − tn−2 )(t − tn−1 )| ≤ = 1. (kn−1 + kn )kn (kn−1 + kn )kn
Hence, the interpolant U ] is superconvergent, uniformly in t, for appropriate γ. Theorem 5.3. Suppose that the time mesh satisfies (4.9), (4.10) and (5.4), and that u satisfies (1.6). Then, with γ ] = 3/σ ] , the postprocessed solution (5.3) satisfies ( ] k γσ , 1 ≤ γ < 3/σ ] , n ] 5 kU − ukJ ≤ (2 + 4 Λ) max kU− − u(tn )k + CT M × 0≤n≤N k3 , γ ≥ 3/σ ] . Proof. Write U ] − u = L(U − u) + (Lu − u) and apply Lemmas 5.1 and 5.2. 6. Spatial discretization. 6.1. The fully discrete DG method. We denote the norm of u in H r (Ω) by kukr , and assume now that A = −∇2 in a bounded, convex or C 2 domain Ω in Rd , subject to homogeneous Dirichlet boundary conditions. Thus, if u ∈ H01 (Ω) and Au ∈ L2 (Ω), then u ∈ H 2 (Ω) and kuk2 ≤ CkAuk (kvk2 := kvkH 2 (Ω) and kvk := kvkL2 (Ω) ). Let Sh ⊆ D(A1/2 ) = H01 (Ω) denote the space of continuous, piecewiselinear functions with respect to a quasi-uniform partition of Ω into triangular or quadrilateral (or tetrahedral etc.) finite elements, with maximum diameter h. Recall that the L2 -projector Ph : L2 (Ω) → Sh and the Ritz projector Rh : H01 (Ω) → Sh are defined by hPh v, W i = hv, W i and A(Rh v, W ) = A(v, W )
for all W ∈ Sh ,
(6.1)
and that the latter has the quasi-optimal approximation property kv − Rh vk + hk∇(v − Rh v)k ≤ Ch2 kvk2
for v ∈ H01 (Ω) ∩ H 2 (Ω).
(6.2)
Let W(Sh ) denote the space of piecewise linear functions U : J → Sh (so U is continuous in space, but may be discontinuous in time). We define the fully discrete DG solution Uh ∈ W(Sh ) by requiring (2.4) to hold for every X ∈ W(Sh ). Equivalently, cf. (2.6), Z tN 0 0 GN (Uh , X) = hUh− , X+ i+ hf (t), X(t)i dt for all X ∈ W(Sh ), (6.3) 0
16
KASSEM MUSTAPHA AND WILLIAM MCLEAN
0 where, for simplicity, we choose Uh− = (Uh )0− = Ph u0 . In view of (2.7), the Galerkin orthogonality property (2.8) now takes the form 0 GN (Uh − u, X) = hPh u0 − u0 , X+ i=0
for all X ∈ W(Sh ).
(6.4)
Similarly, the fully discrete DG solution Zh ∈ W(Sh ) for the dual problem (1.5) is defined by N GN (V, Zh ) = hV−N , Zh+ i for all V ∈ W(Sh ),
N with Zh+ = Ph zT ,
(6.5)
and, since z satisfies (3.3), GN (V, Zh − z) = hV−N , Ph zT − zT i = 0
for all V ∈ W(Sh ).
(6.6)
Theorem 3.3 generalizes as follows. Theorem 6.1. If u and z are the solutions of the initial value problem (1.1) and the dual problem (1.5), and if Uh and Zh are the corresponding fully discrete DG solutions satisfying (6.3) and (6.5), then
N Uh− − u(tN ), zT = GN (u − Π− Rh u, Zh − z) for every zT ∈ H. Proof. By taking V = Uh in (6.5) we find that N N N N hUh− , zT i = hUh− , Ph zT i = hUh− , Zh+ i = GN (Uh , Zh ),
and taking v = u in (3.3) we have hu(tN ), zT i = GN (u, z), so
N Uh− − u(tN ), zT = GN (Uh , Zh ) − GN (u, z) = GN (Uh − u, Zh ) + GN (u, Zh − z) = GN (u, Zh − z), where the final step used (6.4) with X = Zh . Since GN (u, Zh − z) = GN (u − Π− Rh u, Zh − z) + GN (Π− Rh u, Zh − z), the result follows after putting V = Π− Rh u in (6.6). 6.2. Error in the fully discrete DG solution of the dual problem. We modify the splitting (3.6), by writing A−1 (Zh − z) = ζ + Ψ + Φ where ζ = A−1 (Π+ z − z),
Ψ = A−1 Π+ (Ph z − z),
Φ = A−1 (Zh − Π+ Ph z) ∈ W(Sh ).
Theorem 3.8 generalizes as follows. Theorem 6.2. Let z denote the solution of the dual problem (1.5), and let Zh denote the fully discrete DG solution defined by (6.5). Then, for −1 < α < 1, α kA−1 (Zh − z)kJ ≤ C tN+ k 1+α− + h2 `(tN /kN )kzT k. Proof. We already estimated kζk in Lemma 3.5. To estimate Ψ, observe that since A−1 commutes with Π+ and since ARh = Ph A (implying A−1 Ph = Rh A−1 ), Ψ = Π+ (Rh − I)A−1 z.
(6.7)
17
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
Using (2.16), the error bound (6.2) for the Ritz projection, H 2 -regularity for A and Lemma 3.4, we find that kΨkIn ≤ 3k(Rh − I)A−1 zkIn ≤ Ch2 kA−1 z(t)k2 ≤ Ch2 kz(t)k ≤ Ch2 kzT k.
(6.8)
To estimate Φ, let V ∈ W(Sh ) and observe that since A−1 V = A−1 Ph V = Rh A−1 V , GN (V, ζ + Ψ + Φ) = GN (A−1 V, Zh − z) = GN (Rh A−1 V, Zh − z) = 0,
(6.9)
where we used (6.6) with V replaced by Rh A−1 V . From the proof of Lemma 3.6, Z T hV, Bα∗ Aζi dt, (6.10) GN (V, ζ) = − 0
and by (3.2), GN (V, Ψ) =
hV−N , ΨN −i
−
N −1 X
hV−n , [Ψ]n i
N Z X + −hV, Ψ0 i + A(V, Bα∗ Ψ) dt. (6.11)
n=1
n=1
In
Since A(V, Bα∗ Ψ) = A V, Bα∗ A−1 Π+ (Ph − I)z) = hV, (Ph − I)Bα∗ Π+ zi = 0, |GN (V, Ψ)| ≤ kV kJ
kΨN −k
+
N −1 X
n
k[Ψ] k +
n=1
N Z X n=1
kΨ k dt . 0
In
2 By (6.8), kΨN − k ≤ kΨkIN ≤ Ch kzT k. Using (6.7), (2.16) and (6.2) we have Z Z n −1 0 2 k[Ψ] k ≤ k(Rh − I)A z (t)k dt ≤ Ch kA−1 z 0 (t)k2 dt, In
In
so H 2 -regularity for A and the estimate (T − t)kz 0 (t)k ≤ CkzT k from Lemma 3.4 give N −1 X
k[Ψ]n k ≤ Ch2 kzT k
Z
tN −1
(T − t)−1 dt = Ch2 kzT k log(tN /kN ).
0
n=1
Using (6.7), (2.16), (6.2) and Lemma 3.4, we find that Z Z Z 2 Ch2 tn − t kΨ0 (t)k dt ≤ (tn − t)k(Rh − I)A−1 z 0 (t)k dt ≤ kzT k dt, kn In kn In In T − t so N Z X n=1
In
Z kΨ k dt ≤ Ch kzT k 0
2
0
tN −1
−1
(T − t)
dt +
−1 kN
Z
dt ,
IN
and we conclude that |GN (V, Ψ)| ≤ Ch2 `(tN /kN )kV kJ kzT k. Therefore, by (6.9), (6.10) and Lemma 3.7, α |GN (V, Φ)| = |GN (V, ζ)+GN (V, Ψ)| ≤ C tN+ k 1+α− +h2 `(tN /kN )kV kJ kzT k. (6.12) Fix n with 1 ≤ n ≤ N , and define V ∈ W(Sh ) by ( 0, if t ∈ Ij for 1 ≤ j ≤ n − 1, V (t) = Φ(t), if t ∈ Ij for n ≤ j ≤ N .
18
KASSEM MUSTAPHA AND WILLIAM MCLEAN
In view of (2.9), GN (V, Φ) ≥ mate (6.12) gives
1 N 2 2 kΦ− k
2 + 12 kΦn−1 + k +
1 2
PN −1 j=n
k[Φ]j k2 , so the esti-
α+ 1+α− 2 n 2 + h2 `(tN /kN )kΦk(tn−1 ,T ) kzT k kΦn−1 + k + k[Φ] k ≤ C tN k for 1 ≤ n ≤ N − 1, whereas α+ 1+α− −1 2 2 + h2 `(tN /kN )kΦk(tN −1 ,T ) kzT k. kΦN k + kΦN − k ≤ C tN k + n Furthermore, kΦkIn = max kΦn−1 + k, kΦ− k because Φ is piecewise linear in t, and 2 n 2 n 2 kΦn− k ≤ kΦn+ k + k[Φ]n k, implying that kΦk2In ≤ kΦn−1 + k + kΦ+ k + k[Φ] k . By letting n∗ = argmax1≤n≤N kΦkIn , we see that α kΦk2J = kΦk2In∗ ≤ C tN+ k 1+α− + h2 `(tN /kN )kΦkJ kzT k, giving the desired bound for kΦkJ . 6.3. Fully-discrete error . As claimed in the Introduction, we have the following error bound for Uh ; recall that γ ∗ = (2 + α− )/σ whereas γ ] = 3/σ ] . Theorem 6.3. Assume that the solution u of the initial value problem (1.1) satisfies (1.3) and (1.4), and that the time mesh satisfies assumptions (4.9) and (4.10). Then for 0 ≤ n ≤ N the fully discrete DG solution Uh ∈ W(Sh ) satisfies 1+α− +γσ `(tn /kn ), 1 ≤ γ < γ ∗ , k n 2 kUh− − u(tn )k ≤ CT M h + CT M × k 3+2α− `(tn /kn )2 , γ = γ∗, 3+2α− k `(tn /kn ), γ > γ∗. Moreover, if the additional regularity property (1.6) holds, then the fully discrete, postprocessed solution Uh] = LUh satisfies ( ] k γσ , 1 ≤ γ < γ ] , ] n 5 kUh − ukJ ≤ (2 + 4 Λ) max kUh− − u(tn )k + CT M × 0≤n≤N k3 , γ ≥ γ]. Proof. In view of Lemma 4.2, it suffices to show (cf. Theorem 4.1) that n kUh− − u(tn )k ≤ CT k 1+α− + h2 `(tn /kn )(u) + CT M h2 . Put ξ = u − Rh u and η = u − Π− u so that u − Π− Rh u = η + Π− ξ and thus, N by Theorem 6.1, hUh− − u(tN ), zT i = GN (η, Zh − z) + GN (Π− ξ, Zh − z). Using Theorem 6.2 in place of Theorem 3.8, we can show as in the proof of Theorem 4.1 that |GN (η, Zh − z)| ≤ CT kzT k k 1+α− + h2 `(tn /kn ) (u). By (2.5), GN (Π− ξ, Zh ) equals
−
(Π
0 ξ)0+ , Zh+
+
N −1 X
N
− n n X [Π ξ] , Zh+ +
n=1
n=1
Z
(Π− ξ)0 , Zh + A(Bα Π− ξ, Zh ) dt,
In
and, since Bα commutes with the Ritz projector Rh , the definition (6.1) of Rh implies that A(Bα Π− ξ, Zh ) = A(Bα Π− u, Zh ) − A(Rh Bα Π− u, Zh ) = 0. Integrating by parts, n and noting that (Π− ξ)n− = ξ− = ξ(tn ), Z Z
− 0
n n − n−1 n−1
− (Π ξ) , Zh dt = ξ− , Zh− − (Π ξ)+ , Zh+ − Π ξ, Zh0 dt, In
In
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
19
and, because Zh0 is constant on In , the orthogonality property of Π− in (2.10) gives Z Z Z
− n−1 n−1 0 0 n n Π ξ, Zh dt = hξ, Zh i dt = hξ− , Zh− i − hξ+ , Zh+ i − hξ 0 , Zh i dt. In
In
In
Thus, Z In
n−1 n−1 (Π− ξ)0 , Zh dt = ξ+ − (Π− ξ)n−1 + + , Zh+
Z
hξ 0 , Zh i dt
In
and so −1 N
NX
n n X 0 GN (Π− ξ, Zh ) = ξ(0), Zh+ + [ξ] , Zh+ + n=1
Z
hξ 0 , Zh i dt.
In
n=1
Using (3.3) with v = Π− ξ, and noting that [ξ]n = 0, we obtain N
X 0 GN (Π ξ, Zh − z) = ξ(0), Zh+ − ξ(T ), zT +
Z
−
hξ 0 , Zh i dt.
In
n=1
Stability of the fully discrete dual problem, kZh kJ ≤ CkzT k, follows from (2.9), giving Z T 0 GN (Π− ξ, Zh − z) ≤ CkzT k kξ(0)k + kξ(tN )k + kξ k dt 0
Z 2 ≤ Ch kzT k ku0 k2 + ku(T )k2 +
T
ku (t)k2 dt , 0
0 n where we used the error bound (6.2) for the Ritz projector. The error bound for Uh− ] 2 follows using H -regularity of A and the assumption (1.3), and the bound for Uh is proved in the same way as for U ] in Theorem 5.3.
7. Numerical results. We present a series of numerical tests using a model problem in one space dimension, of the form (1.1) with Au = −uxx , Ω = (0, 1), [0, T ] = [0, 1], and homogeneous Dirichlet (absorbing) boundary conditions. These tests reveal faster than expected convergence when α < 0, and that our regularity assumptions are more restrictive than is needed in practice. We apply the fully discrete DG method defined in Section 6.1, employing a time mesh of the form (4.11), for various choices of the mesh grading parameter γ ≥ 1, and a uniform spatial mesh consisting of M subintervals, each of length h = 1/M . We choose M = dN 3/2 e unless indicated otherwise, and always ensure that the spatial error is negligible compared to the error from the time discretization. √ 7.1. Nodal errors. We choose u0 (x) = φ1 (x)/ 2 = sin(πx) for the initial data and f (t, x) = (α + 1) tα sin(πx) for the inhomogeneous term, so that Γ(α + 2) Γ(α + 2) 2 α+1 u(t) = Eα+1 (−π t ) 1− + sin(πx), (7.1) π2 π2 P∞ where the Mittag–Leffler function is given by Eν (t) = p=0 tp /Γ(1 + νp). One easily verifies that the regularity conditions (1.3), (1.4) and (1.6) hold for σ = 1 + α = σ] ,
(7.2)
20
KASSEM MUSTAPHA AND WILLIAM MCLEAN Table 7.1 The left nodal error for Example 1 when α = −0.6, and so γ ∗ = 7/2 = 3.5.
N 20 40 80 160 320 640 ρ
N =M γ=1 6.81e-03 6.02e-03 0.177 5.10e-03 0.238 4.07e-03 0.325 3.06e-03 0.415 2.16e-03 0.498 0.800
γ=3 8.68e-04 2.41e-04 1.848 5.99e-05 2.008 1.38e-05 2.123
M = dN 3/2 e γ = 7/2 5.98e-04 1.50e-04 2.066 3.42e-05 2.171 7.37e-06 2.236
γ=4 4.76e-04 1.14e-04 2.141 2.51e-05 2.216 5.30e-06 2.268
1.600
1.800
1.800
Table 7.2 Left nodal errors when α = −0.25, and so γ ∗ = 7/3.
N 20 40 80 160 ρ
γ=1 5.37e-03 1.83e-03 1.556 5.58e-04 1.712 1.66e-04 1.752 1.500
γ=2 2.57e-04 4.27e-05 2.588 6.78e-06 2.653 1.05e-06 2.685 2.250
γ = 7/3 2.18e-04 3.55e-05 2.712 5.58e-06 2.717 8.62e-07 2.719 2.500
γ=3 2.05e-04 3.34e-05 2.714 5.27e-06 2.711 8.13e-07 2.721 2.500
Table 7.3 Right nodal errors when α = −0.25.
N 20 40 80 160
γ=1 1.31e-01 7.74e-02 0.809 4.42e-02 0.809 2.50e-02 0.821
γ = 7/3 6.69e-03 1.73e-03 2.017 4.91e-04 1.853 1.44e-04 1.789
γ=3 4.94e-03 1.26e-03 2.039 3.20e-04 2.016 8.06e-05 2.008
γ=4 4.82e-03 1.23e-03 2.035 3.14e-04 2.013 7.91e-05 2.006
because Au(t) = π 2 u(t) and A2 u(t) = π 4 u(t). n = limt→t± Uh (t), we use the terminology With Uh± n n left nodal error = max ||Uh− − u(tn )||, 1≤n≤N
n right nodal error = max ||Uh+ − u(tn )||, 1≤n≤N
and expect the former, but not the latter, to exhibit superconvergence. By Theorem 6.3, the left nodal error is of order k ρ (ignoring the logarithmic factors) where ( 1 + α− + γσ, 1 ≤ γ < γ ∗ = (2 + α− )/σ, ρ= (7.3) 3 + 2α− , γ ≥ γ∗, provided h2 k ρ . Tables 7.1 and 7.2 show the left nodal errors when α = −0.6 and −0.25, respectively, for various choices of N and γ. We observe that the convergence rates are faster than predicted by our theory, for σ given by (7.2), and appear to be O(k 3+α ) when γ ≥ γ ∗ . Table 7.3 shows the right nodal errors when α = −0.25, which never exhibit better than O(k 2 ) convergence, and thus are not superconvergent. Tables 7.4 and 7.5 show the left and right nodal errors, respectively, when α = +0.25, and Table 7.6 shows the left nodal errors when α = +0.75. For these positive
21
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD Table 7.4 Left nodal errors when α = +0.25, and so γ ∗ = 8/5 = 1.6.
N 20 40 80 160 ρ
γ=1 4.98e-05 1.22e-05 2.103 2.57e-06 2.288 4.91e-07 2.408 2.250
γ = 1.3 4.33e-05 5.20e-06 3.170 6.23e-07 3.116 7.54e-08 3.074 2.625
γ = 8/5 4.45e-05 5.38e-06 3.157 6.52e-07 3.102 7.94e-08 3.065 3.000
γ = 1.9 5.11e-05 6.22e-06 3.146 7.60e-07 3.087 9.34e-08 3.053 3.000
Table 7.5 Right nodal errors when α = +0.25.
N 20 40 80 160
γ=1 1.00e-02 5.22e-03 0.973 2.39e-03 1.149 1.04e-03 1.211
γ = 1.3 4.17e-03 1.46e-03 1.572 4.84e-04 1.619 1.58e-04 1.628
γ = 8/5 3.04e-03 7.67e-04 2.059 1.93e-04 2.027 4.84e-05 2.015
γ = 1.9 3.27e-03 8.33e-04 2.044 2.09e-04 2.027 5.25e-05 2.015
Table 7.6 Left nodal errors when α = +0.75, and so γ ∗ = 8/7 ≈ 1.143.
N 20 40 80 160 ρ
γ=1 5.85e-05 7.51e-06 3.067 9.46e-07 3.042 1.19e-07 3.018 2.750
γ = 1.1 5.48e-05 6.99e-06 3.080 8.77e-07 3.047 1.10e-07 2.995 2.925
γ = 8/7 5.42e-05 6.90e-06 3.081 8.65e-07 3.048
γ = 1.3 5.34e-05 7.00e-06 3.036 8.95e-07 3.022
3.000
3.000
values of α, the left nodal errors appear to be O(k 3 ) when γ is sufficiently large, whereas the right nodal errors are again at best O(k 2 ). 7.2. Uniform error after post-processing. We introduce a finer mesh G N,m = { tj−1 + `kj /m : j = 1, 2, . . . , N and ` = 0, 1, . . . , m },
(7.4)
and define the discrete maximum norm kvkJ,m = maxt∈G N,m kv(t)k, so that, for sufficiently large values of m, kUh] − ukJ,m approximates the uniform error kUh] − ukJ , ] which, by Theorem 6.3, is of order k ρ where ( min(ρ, γσ ] ), 1 ≤ γ < γ ] = 3/σ ] , ] ρ = min(ρ, 3), γ ≥ γ]. Tables 7.7 and 7.8 show kUh] − ukJ,12 when α = −0.25 and α = +0.25, respectively. Comparing these results with those of Tables 7.2–7.5 confirms that the full rate of nodal superconvergence is reflected in the uniform error after postprocessing, but only with a stronger mesh grading, since γ ] > γ ∗ . 7.3. A less regular solution. Now let u0 (x) = x(1 − x) and f ≡ 0. Thus, Au0 = 2 does not vanish at 0 and 1, so u0 ∈ / D(A2 ) [19, Equation (3.2)]. Separation of variables yields a series representation u(x, t) = 8
∞ X n=0
ωn−3 sin(ωn x)E1+α (−ωn2 t1+α )
with ωn = (2n + 1)π,
(7.5)
22
KASSEM MUSTAPHA AND WILLIAM MCLEAN Table 7.7 Uniform errors after postprocessing when α = −0.25, so γ ∗ = 7/3 and γ ] = 4.
N 20 40 80 160 320 ρ]
γ=2 8.81e-03 2.95e-03 1.635 1.02e-03 1.559 3.58e-04 1.525 1.26e-04 1.505 1.500
γ = 7/3 4.98e-03 1.30e-03 2.006 3.67e-04 1.859 1.07e-04 1.791 3.17e-05 1.755 1.750
γ=3 2.76e-03 5.48e-04 2.418 1.14e-04 2.308 2.39e-05 2.274 5.01e-06 2.253 2.250
γ=4 1.75e-03 2.14e-04 3.140 2.71e-05 3.032 3.51e-06 2.976 4.60e-07 2.931 2.500
Table 7.8 Uniform errors after postprocessing when α = +0.25, and so γ ∗ = 8/5 = 1.6 and γ ] = 12/5 = 2.4.
N 20 40 80 160 ρ]
γ = 1.3 3.11e-03 1.07e-03 1.590 3.55e-04 1.625 1.16e-04 1.630 1.625
γ = 8/5 1.08e-03 3.17e-04 1.827 8.28e-05 1.972 2.09e-05 2.002 2.000
γ = 1.9 6.49e-04 1.35e-04 2.353 2.63e-05 2.398 5.08e-06 2.393 2.375
γ = 12/5 6.43e-04 7.89e-05 3.135 9.63e-06 3.089 1.18e-06 3.052 3.000
and it is not difficult to verify directly that (1.3) holds. In fact, differentiating (7.5), ∂tj Au(x, t) = ∂tj uxx = −8
∞ X
ωn−1 sin(ωn x)
n=0
dj E1+α (−ωn2 t1+α ), dtj
and thus by Parseval’s identity, k∂tj Au(t)k2 = 32
∞ X
ωn−2
n=0
2 dj 2 1+α E (−ω t ) . 1+α n dtj
P∞ Since |E1+α (−ωn2 t1+α )| ≤ C for all t > 0 it follows that kAu(t)k2 ≤ C n=0 ωn−2 < ∞, and the Mittag–Leffler function satisfies [19, Theorem 4.2] j d −(1+α)µ−j −2µ 2 1+α ωn for j ∈ {1, 2, 3, . . .} and |µ| ≤ 1, (7.6) E (−ω t ) n dtj 1+α ≤ Ct so j
t
2 k∂tj Au(t)k
−2µ(1+α)
≤ Ct
∞ X
ωn−4µ−2
n=0
C t−µ(1+α) ≤ 1 + 4µ
2 for − 14 < µ ≤ 1.
Thus, the regularity conditions (1.3) hold for σ < 41 (1+α) and M = C(1+α−4σ)−1/2 . Similarly, ∂tj A2 u(x, t) = ∂tj uxxxx = 8
∞ X
ωn sin(ωn x)
n=0
and
k∂tj A2 u(t)k2
= 32
P∞
2 n=0 ωn
dj E1+α (−ωn2 t1+α ) dtj
2
dj 2 1+α ) dtj E1+α (−ωn t
, so
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
23
Table 7.9 Left nodal errors for the less regular solution (7.5).
N 20 40 80 160 240 ρ
γ=1 1.84e-03 1.68e-03 0.129 1.44e-03 0.223 1.15e-03 0.320 9.84e-04 0.393 0.500
N 20 80 160 ρ
γ=1 2.01e-03 3.90e-04 1.143 2.21e-04 0.821 0.875
N 20 40 80 160 ρ
γ=1 2.10e-04 6.77e-05 1.632 2.19e-05 1.636 7.11e-06 1.625 1.025
N 20 40 80 160
γ=1 7.37e-05 1.83e-05 2.011 4.56e-06 2.003 1.14e-06 2.001
j
t
2 k∂tj A2 u(t)k
α = −0.6, γ ∗ = 14 γ=3 γ=5 2.47e-04 1.09e-04 6.99e-05 1.818 2.54e-05 2.103 2.26e-05 1.628 5.52e-06 2.202 8.06e-06 1.489 1.16e-06 2.250 4.38e-06 1.500 4.60e-07 2.282 0.700 0.900 ∗ α = −0.3, γ = 68/7 ≈ 9.71 γ=2 γ=3 1.08e-04 6.39e-05 9.33e-06 1.758 1.81e-06 2.596 2.77e-06 1.753 2.92e-07 2.632 1.050 1.225 α = +0.3, γ ∗ = 80 γ = 1.5 γ = 1.75 2.08e-05 1.21e-05 3.61e-06 2.527 1.61e-06 2.904 6.43e-07 2.486 2.13e-07 2.917 1.17e-07 2.461 2.80e-08 2.930 1.038 1.044 α = +0.6 γ = 1.25 γ = 1.5 1.88e-05 1.30e-05 3.12e-06 2.590 1.87e-06 2.801 5.31e-07 2.554 2.60e-07 2.845 9.21e-08 2.528 3.64e-08 2.840
−2(1+α)µ
≤ Ct
∞ X n=0
ωn2−4µ ≤
C(t−(1+α)µ )2 4µ − 3
γ=6 1.05e-04 2.44e-05 2.111 5.29e-06 2.203 1.11e-06 2.247 2.30e-07 2.278 1.000 γ = 3.25 6.39e-05 1.82e-06 2.595 2.94e-07 2.632 1.269 γ=2 1.23e-05 1.57e-06 2.966 1.99e-07 2.983 2.53e-08 2.972 1.050 γ=2 2.19e-05 3.14e-06 2.803 4.22e-07 2.894 5.51e-08 2.938
for
3 4
< µ ≤ 1.
Thus, (1.4) holds for σ < 1 + α− − 34 (1 + α), but σ > 0 only if −1 < α < 31 . Table 7.9 shows the left nodal errors for α = −0.6, −0.3, +0.3 and +0.6. To compute the error, we evaluated u by truncating the Fourier series (7.5) after 160 terms. In the first three cases, α < 13 and we show the values of ρ given by (7.3) with σ = 1 + α− − 43 (1 + α). Evidently, the theoretical convergence rates are even more pessimistic than in our previous example (7.1), owing to the smaller value of σ, and the values of γ ∗ become unrealistically large. However, in practice, we observe that the errors are superconvergent of order k 3+α− for moderate γ, even in the case α = +0.6 when the regularity assumption (1.4) fails completely. 8. Concluding remarks. We have analysed a piecewise-linear DG method for the time discretization of (1.1) — a fractional diffusion (−1 < α < 0) or wave (0 < α < 1) equation — and proved superconvergence at the nodes, generalizing a known result for the classical heat equation. Numerical experiments indicate that our theoretical error bounds are sharp if α > 0, but not if α < 0. For generic regular
24
KASSEM MUSTAPHA AND WILLIAM MCLEAN
data u0 and f , derivatives of the exact solution are singular as t → 0, but nevertheless by employing non-uniform time steps we achieve a high convergence rate of O(k 3+α− ). After postprocessing the solution, the same high accuracy is achieved for all t, not just at the nodes. We have also proved that the additional error arising from a spatial discretization by continuous piecewise-linear finite elements is O(h2 ). REFERENCES [1] C.-M. Chen, F. Liu, V. Anh, and I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput., 32 (2010), pp. 1740–1760. , Numerical methods for solving a two-dimensional variable-order anomalous subdiffu[2] sion equation, Math. Comp., 81 (2012), pp. 345–366. [3] E. Cuesta, C. Lubich, and C. Palencia, Convolution quadrature time discretization of fractional diffusion-wave equations, Math. Comp., 75 (2006), pp. 673–696. [4] E. Cuesta and C. Palencia, A fractional trapezoidal rule for integro-differential equations of fractional order in Banach spaces, Appl. Numer. Math., 45 (2003), pp. 139–159. , A numerical method for an integro-differential equation with memory in Banach spaces: [5] qualitative properties, SIAM J. Numer. Anal., 41 (2003), pp. 1232–1241. [6] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys., 228 (2009), pp. 7792–7804. [7] M. Cui, Compact alternating direction implicit method for two-dimensional time fractional diffusion equation, J. Comput. Phys., 231 (2012), pp. 2621–2633. [8] , Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algor., (Published online: 2012). ´e, Time discretization of parabolic problems by the [9] K. Eriksson, C. Johnson, and V. Thome discontinuous Galerkin method, M2AN Math. Model. Numer. Anal., 19 (1985), pp. 611– 643. [10] A. Hanyga, Wave propagation in media with singular memory, Math. Comput. Modelling., 34 (2001), pp. 1399–1421. [11] B. Jin, R. Lazarov, and Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations. Preprint, arXiv:1204.38884v1. [12] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and Applications of Fractional Differential Equations, vol. 204 of North-Holland Mathematics Studies, North–Holland, Amsterdam, 2006. [13] T. A. M. Langlands and B. I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys., 205 (2005), pp. 719–736. [14] F. Liu, C. Yang, and K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math., 231 (2009), pp. 160–176. ´ pez-Ferna ´ ndez and C. Palencia, On the numerical inversion of the Laplace transform [15] M. Lo of certain holomorphic functions, Appl. Numer. Math., 51 (2004), pp. 289–303. ´ pez-Fernandez, C. Palencia, and A. Scha ¨ dle, A spectral order method for inverting [16] M. Lo sectorial Laplace transforms, SIAM J. Numer. Anal., 44 (2006), pp. 1332–1350. [17] F. Mainardi and P. Paradisi, Fractional diffusive waves, J. Comput. Acoustics, 9 (2001), pp. 1417–1436. [18] A. M. Mathai, R. K. Saxena, and H. J. Haubold, The H-Function: Theory and Applications, Springer, 2010. [19] W. McLean, Regularity of solutions to a time-fractional diffusion equation, ANZIAM J., 52 (2010), pp. 123–138. [20] , Fast summation by interval clustering for an evolution equation with memory. Preprint, arXiv:1203.4032v1, 2012. [21] W. McLean and K. Mustapha, A second-order accurate numerical method for a fractional wave equation, Numer. Math., 105 (2007), pp. 481–510. [22] , Convergence analysis of a discontinuous Galerkin method for a sub-diffusion equation, Numer. Algor., 52 (2009), pp. 69–88. ´e, Numerical solution of an evolution equation with a positive-type [23] W. McLean and V. Thome memory term, ANZIAM J., 35 (1993), pp. 23–70. ´e, Numerical solution via Laplace transforms of a fractional order [24] W. McLean and V. Thome
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
25
evolution equation, J. Integral Equations Appl., 22 (2010), pp. 57–94. ´e, and L. B. Wahlbin, Discretization with variable time steps of an [25] W. McLean, V. Thome evolution equation with a positive-type memory term, J. Comput. Appl. Math., 69 (1996), pp. 49–69. [26] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Reports, 339 (2000), pp. 1–77. , The restaurant at the end of the random walk: Recent developments in the description [27] of anomalous transport by fractional dynamics, J. Phys. A, 37 (2004), pp. R161–R208. [28] K. Mustapha, An implicit finite difference time-stepping method for a sub-diffusion equation, with spatial discretization by finite elements, IMA Journal Numer. Anal., 31 (2011), pp. 719–739. [29] K. Mustapha and J. AlMuttawa, A finite difference method for an anomalous sub-diffusion equation, theory and applications, Numer. Algor., (2012). [30] K. Mustapha and W. McLean, Discontinuous Galerkin method for an evolution equation with a memory term of positive type, Math. Comp., 78 (2009), pp. 1975–1995. , Piecewise-linear, discontinuous Galerkin method for a fractional diffusion equation, [31] Numer. Algor., 56 (2011), pp. 159–184. , Uniform convergence for a discontinuous Galerkin, time stepping method applied to a [32] fractional diffusion equation, IMA J. Numer. Anal., 32 (2012), pp. 906–925. [33] Y. nan Zhang and Z. zhong Sun, Alternating direction implicit schemes for the twodimensional fractional sub-diffusion equation, J. Comput. Phys., 230 (2011), pp. 8713– 8728. [34] I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, 1999. ¨ ss, Evolutionary Integral Equations and Applications, vol. 87 of Monographs in Mathe[35] J. Pru matics, Birkh¨ auser, Basel, 1993. [36] J. M. Sanz-Serna, A numerical method for a partial integro-differential equation, SIAM J. Numer. Anal., 25 (1988), pp. 319–327. ¨ dle, M. Lo ´ pez-Ferna ´ ndez, and C. Lubich, Fast and oblivious convolution quadra[37] A. Scha ture, SIAM J. Sci. Comput., 28 (2006), pp. 421–438. ´ ndez, and R. J. Cherry, Anoma[38] P. R. Smith, I. E. G. Morrison, K. M. Wilson, N. Ferna lous diffusion of major histocompatability complex class I molecules on hela cells determined by single particle tracking, Biophys. J., 76 (1999), pp. 3331–3344. [39] I. Sokolov and J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos, 15 (2005), p. 026103. [40] V. E. Tarasov, Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media (Nonlinear Physical Science), Springer, 2010. [41] H. Wang and K. Wang, An O(N log2 N ) alternating-direction finite difference method for two-dimensional fractional diffusion equations, J. Comput. Phys., 230 (2011), pp. 7830– 7839. [42] S. B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys., 216 (2006), pp. 264–274. [43] S. B. Yuste and L. Acedo, An explicit finite difference method and a new von Neumanntype stability analysis for fractional diffusion equations, SIAM J. Numer. Anal., 42 (2005), pp. 1862–1874. [44] P. Zhuang, F. Liu, V. Anh, and I. Turner, New solution and analytical techniques of the implicit numerical methods for the anomalous sub-diffusion equation, SIAM J. Numer. Anal., 46 (2008), pp. 1079–1095. [45] , Stability and convergence of an implicit numerical method for the nonlinear fractional reaction-subdiffusion process, IMA J. Appl. Math., 74 (2009), pp. 645–667.
1. Introduction. In previous work [22, 30, 31, 32], we have studied discontinuous Galerkin (DG) methods for the time discretization of the abstract intial value problem u0 + Bα Au = f (t)
for 0 < t < T ,
with u(0) = u0 ,
(1.1)
where u0 = ∂u/∂t and Bα = ∂t−α ; more precisely, letting ωµ (t) = tµ−1 /Γ(µ) for µ > 0, the function Bα v is either a (Riemann–Liouville) fractional order derivative in time, ∂ Bα v(t) = ∂t
t
Z
ω1+α (t − s)v(s) ds if −1 < α < 0,
(1.2)
0
or a fractional order integral in time, Z Bα v(t) =
t
ωα (t − s)v(s) ds if 0 < α < 1. 0
In Section 2 we set out technical assumptions on the operator A, but for the present discussion we simply take Au = −∇2 u on a spatial domain Ω, and impose homogeneous Dirichlet boundary conditions on u. Problems of the form (1.1) arise in a variety of physical, biological and chemical applications [12, 18, 26, 27, 34, 38, 39, 40]. The case −1 < α < 0 describes slow or anomalous sub-diffusion and occurs, for example, in models of fractured or porous media, where the particle flux depends on the entire history of the density gradient ∇u. The case 0 < α < 1 describes wave propagation in viscoelastic materials [10, 17, 35]. ∗ We
thank KFUPM for supporting this research as part of the project SB101020. Department of Mathematics and Statistics, KFUPM, Dhahran 31261, Saudi Arabia ([email protected]). ‡ School of Mathematics and Statistics, The University of New South Wales, Sydney 2052, Australia ([email protected]). †
1
2
KASSEM MUSTAPHA AND WILLIAM MCLEAN
In the limit as α → 0, the evolution equation in (1.1) becomes u0 + Au = f , which is just the classical heat equation, and Eriksson et al. [9] studied the convergence of the DG solution U (t) ≈ u(t) in this case. For a maximum time step k, and using discontinuous piecewise polynomials of degree at most q − 1 in t, with no spatial discretization, they proved for 0 ≤ t ≤ T an optimal convergence rate Z t kf (q) (s)k ds , kU (t) − u(t)k ≤ Ck q ku0 kq + ku(q) (0)k + kf (q−1) (0)k + 0
where kvk is the norm of v in L2 (Ω) and where kvkq = kAq/2 vk for v in D(Aq/2 ), the domain of the unbounded operator Aq/2 . In addition, they proved that the DG solution is superconvergent at the nth time level tn , satisfying an error bound Z tn (q) − 2q−1 (q) kf (s)kq−1 ds , kU (tn ) − u(tn )k ≤ Ck ku0 k2q−1 + ku (0)kq−1 + 0
U (t) denotes the limit from the left. Ericksson et al. were also where U (t− n ) = limt→t− n able to prove that a convergence rate faster then O(k q ) holds under less restrictive spatial regularity requirements on the solution u. Our aim is to establish superconvergence results for the fractional-order problem (1.1), restricting our attention to the piecewise-linear DG method (q = 2). We believe that this time-stepping scheme is the first to achieve better than second-order accuracy. As well as nodal superconvergence of the DG solution we show that a postprocessed solution is superconvergent uniformly in t. Many authors have studied numerical methods for (1.1). In the case 0 < α < 1, Sanz-Serna [36] proposed a convolution quadrature scheme, and subsequently Cuesta, Lubich and Palencia [4, 5, 3] developed this approach to obtain an O(k 2 ) method as well as a fast implementation [37]. McLean and Thom´ee [23] combined finite differences and quadrature in time, with finite elements in space. In the case −1 < α < 0, Langlands and Henry [13] introduced an implicit Euler scheme involving the Gr¨ unwald–Letnikov fractional derivative and spatial finite differences with step size h, and observed O(k 1/2 + h2 ) convergence in the case α = −1/2. Yuste and Acedo [43] treated an explicit Euler scheme and showed O(k + h2 ) convergence. Zhuang, Liu, Anh, Turner et al. [2, 14, 44, 45] developed another class of O(k + h2 ) finite difference methods, and Yuste [42] presented an O(k 2 + h2 ) method. Cui [6] and Chen et al. [1] studied O(k + h4 ) schemes, and Cui [7, 8] analysed an O(k min(1−α,2+α) + h4 ) ADI scheme on a rectangular spatial domain; see also Wang and Wang [41] and Zhang and Sun [33]. For another type of finite difference scheme [28, 29], the error is O(k 2+α + h2 ), and recently Jin et al. [11] proved optimal error bounds for two semidiscrete finite element methods. Some of these works employ an alternative formulation of (1.1) using the Caputo fractional derivative. In practice, the higher order derivatives of u are typically singular [19, 21] as t → 0, so formally high order methods [1, 2, 6, 7, 8, 14, 33, 41, 43, 44, 45] can fail to achieve fast convergence. We have analysed several methods that allow for the singular behaviour of u by employing non-uniform time steps [21, 25, 28, 29, 32]. Another approach, that yields a parallel in time algorithm with spectral accuracy even for problems with low regularity, is to approximate u via the Laplace inversion formula [15, 16, 24]. To minimise the need for handling separately the cases α < 0 and α > 0, it is convenient to write α+ = max(α, 0) ≥ 0 and α− = min(α, 0) ≤ 0 for the positive and
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
3
negative parts of α, respectively. To prove our superconvergent nodal error bounds, we assume that there exist positive constants M and σ such that kAu0 k + kAu(t)k ≤ M
and kAu0 (t)k + tkAu00 (t)k ≤ M tσ−1 ,
(1.3)
as well as tkA2 u0 (t)k + t2 kA2 u00 (t)k ≤ M tσ−α− −1 ,
(1.4)
for 0 < t ≤ T . For instance [19, 21], if f ≡ 0 and u0 ∈ D(A2 ), then (1.3) and (1.4) hold with M = CkAu0 k and σ = 1 + α− . Section 2 sets out our notation and assumptions, and recalls some tools and results from earlier work [31]. In Section 3, we introduce the homogeneous dual problem, −z 0 + Bα∗ Az = 0
for 0 < t < T ,
with z(T ) = zT ,
(1.5)
for a given terminal value zT , and represent the nodal error U (t− n ) − u(tn ) in terms of z(t) and its DG approximation Z(t). We allow a class of non-uniform meshes, specified in Section 4, where we prove in Theorem 4.3 that the nodal error is O(k 3+2α− ). Our method of analysis allows us to handle the two cases −1 < α < 0 and 0 < α < 1 together, but the former presents additional technical difficulties in some places. In an earlier paper [30, Theorem 4.1], we estimated the nodal error for the case 0 < α < 1 in a different way that yields a bound of order k 2+α . (Although we claimed O(k 3 ) convergence, the first line of [30, Corollary 4.2] contains an error.) In Section 5 we construct, via a simple interpolation scheme, a postprocessed solution U ] whose error is O(k 3+2α− ) for all t, not just at the nodal values, assuming that ku0 (t)k + t2 ku000 (t)k ≤ M tσ
]
−1
for 0 < t ≤ T .
(1.6)
For instance [19], if u0 ∈ D(A) and f ≡ 0, then (1.6) holds with σ ] = 1 + α and M = kAu0 k. Section 6 introduces a fully discrete scheme by applying a continuous piecewise-linear, finite element method for the spatial discretization. Thus, the fully discrete solution is continuous in space but discontinuous in time. We show that the error bound is as for the semidiscrete method but with an extra term of order h2 . Finally, we present some numerical examples in Section 7, which indicate that our error bounds are pessimistic, at least in some cases. We observe that the nodal error from the time discretization is O(k 3+α− ), which is better than our theoretical estimate by a factor k α− . The same is true for the postprocessed solution, uniformly in t. 2. Preliminaries. 2.1. Assumptions on the spatial operator. We assume as in earlier work [9, 31] that the self-adjoint linear operator A has a complete eigensystem in a real Hilbert space H, say Aφj = λj φj for j = 1, 2, 3, . . . , and that A is strictly positive-definite with the eigenvalues ordered so that 0 < λ1 ≤ λ2 ≤ · · · . (Strict positive definiteness is not essential, but allowing λ1 = 0 would result in some technical complications [30, Section 5] that we prefer to avoid here.) We denote p the inner product of u and v in H by hu, vi and the corresponding norm by kuk = hu, ui. Associated with the linear operator A is a bilinear form, denoted by the same symbol: A(u, v) =
∞ X m=1
λm hu, φm ihφm , vi for u, v ∈ D(A1/2 ).
4
KASSEM MUSTAPHA AND WILLIAM MCLEAN
These assumptions hold, in particular, if A = −∇2 subject to homogenous Dirichlet boundary conditions on R a bounded domain Ω, because A has a compact inverse on H = L2 (Ω) and A(u, v) = Ω ∇u · ∇v dx. 2.2. The discontinuous Galerkin time discretization. Fixing a time interval [0, T ], we introduce a mesh for the time discretization, 0 = t0 < t1 < · · · < tN = T,
(2.1)
with kn = tn − tn−1 and In = (tn−1 , tn ) for 1 ≤ n ≤ N , and a maximum time step k = max1≤n≤N kn . Let Pr denote the space Sn of polynomials of degree at most r with coefficients in D(A1/2 ), and let Jn = j=1 Ij = [0, tn ] \ {t0 , t1 , . . . , tn }, with J = JN . Our trial space W consists of the piecewise-linear functions U : J → D(A1/2 ) with U |In ∈ P1 for 1 ≤ n ≤ N . We treat U as undefined at each time level tn , and write n U− = U (t− n ),
n U+ = U (t+ n ),
n n [U ]n = U+ − U− .
(2.2)
For r ∈ {0, 1, 2, . . .} we let C r (J, H) denote the space of functions v : J → H such that the restriction v|In extends to an r-times continuously differentiable function on the closed interval [tn−1 , tn ], for 1 ≤ n ≤ N . In other words, v is a piecewise C r function with respect to the time levels tn . If v ∈ C 1 (J, H), then its fractional derivative (1.2) admits the representation [31] Bα v(t) =
0 ω1+α (t)v+
+
n−1 X
j
ω1+α (t − tj )[v] +
j=1
n Z X j=1
min(tj ,t)
ω1+α (t − s)v 0 (s) ds (2.3)
tj−1
for t ∈ In and −1 < α < 0. Thus, Bα v(t) is left-continuous at t = tn−1 but has a n−1 weak singularity (t − tn−1 )α as t → t+ 6= 0. However, if 0 < α < 1 then n−1 if [v] Bα v(t) is continuous for 0 ≤ t ≤ T . For −1 < α < 1, the piecewise-linear DG time stepping procedure determines 0 U ∈ W by setting U− = u0 and requiring [30, 31] n−1 n−1 i , X+ hU+
Z +
0 hU (t), X(t)i + A Bα U (t), X(t) dt
In
=
n−1 n−1 i , X+ hU−
Z hf (t), X(t)i dt, (2.4)
+ In
for 1 ≤ n ≤ N and for every test function X ∈ P1 . The nonlocal nature of the operator Bα means that at each time step we must compute a sum involving all previous times levels, but this sum can be evaluated via a fast algorithm [20]. 2.3. Galerkin orthogonality and stability. For v ∈ C 1 J, D(A1/2 ) and w ∈ C J, D(A1/2 ) , we define the global bilinear form GN (v, w) =
0 0 hv+ , w+ i
+
N −1 X
h[v]
n
n , w+ i
n=1
N Z X 0 + hv , wi + A(Bα v, w) dt. n=1
(2.5)
In
Summing the equations (2.4) gives 0 0 GN (U, X) = hU− , X+ i+
Z
tN
hf (t), X(t)i dt for all X ∈ W, 0
(2.6)
5
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
and conversely, (2.6) implies that U satisfies (2.4) for 1 ≤ n ≤ N . Since [u]n = 0, Z tN 0 GN (u, X) = hu0 , X+ i + hf (t), X(t)i dt, (2.7) 0 0 and thus, assuming U− = u0 , the error has the Galerkin orthogonality property
GN (U − u, X) = 0
for all X ∈ W.
(2.8)
The DG method is unconditionally stable. Indeed, with the notation kU kI = sup kU (t)k
for any I ⊆ [0, T ],
t∈I
the following estimate holds. 0 Theorem 2.1. Given U− ∈ H and f ∈ L1 (0, T ); H , there exists a unique U ∈ W satisfying (2.4) for n = 1, 2, . . . , N . Furthermore, U (t) ∈ D(A) (t > 0), and Z tn 0 0 2 kU kJn ≤ 8 hU− , U+ i + hf (t), U (t)i dt for 1 ≤ n ≤ N . 0
Proof. Since hV 0 , V i = 12 (d/dt)kV k2 and 2GN (V, V ) ≥ kV+0 k2 + kV−N k2 +
RT 0
N −1 X
A(Bα V, V ) dt ≥ 0 we find that k[V ]n k2
for all V ∈ W,
(2.9)
n=1
implying the stated estimate; for details, see [30, Theorem 2.1], [31, Theorem 1]. 2.4. A discontinuous quasi-interpolant. The conditions Z − Π− v(t− ) = v(t ) and v(t) − Π− v(t)] dt = 0 n n
(2.10)
In
determine a unique projection operator Π− : C(J, H) → W; explicitly, Π− v(t) := v(t− n)+
¯n v(t− n)−v (t − tn ) kn /2
for t ∈ In ,
R where v¯n = kn−1 In v(t) dt denotes the mean value of v over In . The interpolation error admits the integral representations [30, Equation (3.8)] Z tn Z tn − t Π− v(t) − v(t) = v 0 (s) ds − 2 (s − tn−1 )v 0 (s) ds 2 k t In n (2.11) Z tn Z t − t n 00 2 00 = (t − s)v (s) ds + (s − tn−1 ) v (s) ds, for t ∈ In . kn2 t In Likewise, the conditions + Π+ v(t+ n−1 ) = v(tn−1 )
Z and
v(t) − Π+ v(t)] dt = 0,
In
determine a unique projector Π+ : C(J, H) → W, with Π+ v(t) := v(t+ n−1 ) +
v¯n − v(t+ n−1 ) (t − tn−1 ) kn /2
for t ∈ In ,
(2.12)
6
KASSEM MUSTAPHA AND WILLIAM MCLEAN
and t
Z t − tn−1 v (s) ds + 2 Π v(t) − v(t) = − (tn − s)v 0 (s) ds kn2 tn−1 In Z Z t t − tn−1 = (s − tn−1 )v 00 (s) ds + (tn − s)2 v 00 (s) ds. 2 k tn−1 In n Z
+
0
(2.13)
Thus, for t ∈ In (Π+ v − v)0 (t) = −v 0 (t) +
2 kn2
Z
(tn − s)v 0 (s) ds
(2.14)
In
and short calculations lead to the error bounds Z kv (r) (t)k dt ≤ (4−r)knr kv (r) kIn kΠ± v −vkIn ≤ (4−r)knr−1
for r ∈ {1, 2}, (2.15)
In
and the stability estimates Z
± 0
(Π ) v ≤ 2 In kn
kΠ± vkIn ≤ 3kvkIn ,
kv 0 (t)k dt,
± n
[Π v] ≤
Z
kv 0 (t)k dt.
In
In
(2.16) 3. Dual problem. 3.1. Properties of the adjoint operator. The adjoint operator appearing in the dual problem (1.5) should satisfy, for appropriate u and v, the identity Z
T
Z
T
hv, Bα wi dt = 0
hBα∗ v, wi dt,
(3.1)
0
and the next lemma establishes an explicit representation of Bα∗ . Lemma 3.1. The identity (3.1) holds in the following cases. 1. If −1 < α < 0 and v, w ∈ C 1 (J, H), with Bα∗ w(t) = −
∂ ∂t
T
Z
ω1+α (s − t)w(s) ds
for t ∈ J.
t
2. If 0 < α < 1 and v, w ∈ C(J, H), with Bα∗ w(t)
Z
T
ωα (s − t)w(s) ds
=
for 0 ≤ t ≤ T .
t
Proof. In case 1, we see from the representation (2.3) that Z
T
hBα v, wi dt = 0
where, letting B n,j =
N Z X n=1
R In
hBα v, wi dt = S1 + S2 + S3 ,
In
ω1+α (t − tj )w(t) dt and
Z
Z
min(tj ,t)
Dnj = In
tj−1
ω1+α (t − s)hv 0 (s), w(t)i ds dt,
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
7
we define S1 =
N X
0 n,0 v+ , B ,
S2 =
n=1
N n−1 X X
[v]j , B n,j ,
S3 =
n=2 j=1
N X n X
Dnj .
n=1 j=1
By reversing the order of integration, integrating by parts and then interchanging variables, we find that for 1 ≤ j ≤ n − 1, Z Z ∂ j j−1 ω1+α (s − t)w(s) ds dt, Dnj = hv− , B n,j i − hv+ , B n,j−1 i − v(t), ∂t In Ij whereas n−1 Dnn = −hv+ , B n,n−1 i −
Z
∂ ∂t
v(t),
In
tn
Z
ω1+α (s − t)w(s) ds dt.
t
Thus, after interchanging the order of summation for the double integrals, N Z N n−1 N X j X X X j−1 n−1 n,j n,j−1 n,n−1 hv, Bα∗ wi dt, hv− , B i − hv+ , B i − i+ S3 = − hv+ , B n=2 j=1
n=1
j=1
Ij
that is, S3 =
N n−1 X X
j 0 hv− , B n,j i − hv+ , B 1,0 i −
n=2 j=1
Z N n−1 X X j hv+ , B n,j i −
hv, Bα∗ wi dt
0
n=2 j=0
Z
T
T
hv, Bα∗ wi dt,
= −S1 − S2 − 0
so (3.1) holds. In the case 0 < α < 1, we simply reverse the order of integration. The adjoint operator admits a representation analogous to (2.3). Lemma 3.2. If −1 < α < 0, then Bα∗ w(t) equals ω1+α (tn −
N t)w−
−
N −1 X
j
ω1+α (tj − t)[w] −
N Z X j=n
j=n
tj
ω1+α (s − t)w0 (s) ds
max(tj−1 ,t)
Bα∗ w(t) is n t− n if [w]
1
for w ∈ C (J, H) and t ∈ In . Thus, right-continuous at t = tn but possesses a weak singularity (tn − t)α as t → 6= 0. Proof. If t ∈ In and j ≥ n + 1, then, integrating by parts, Z Z j j−1 ω1+α (s−t)w(s) ds = ω2+α (tj −t)w− −ω2+α (tj−1 −t)w+ − ω2+α (s−t)w0 (s) ds Ij
Ij
and Z
tn
ω1+α (s − t)w(s) ds = ω2+α (tn −
n t)w−
Z
tn
−
t
ω2+α (s − t)w0 (s) ds.
t
Differentiating these expressions with respect to t, we see from part 1 of Lemma 3.1 that Bα∗ w(t) equals N X j=n
j ω1+α (tj − t)w−
−
N X j=n+1
j−1 ω1+α (tj−1 − t)w+
−
N Z X j=n
tj
max(tj−1 ,t)
and the result follows after shifting the index in the second sum.
ω1+α (s − t)w0 (s) ds,
8
KASSEM MUSTAPHA AND WILLIAM MCLEAN
3.2. Representation of the nodal error. Integration by parts in (2.5), together with the identity (3.1), shows that for all v, w ∈ C 1 (J, H), GN (v, w) =
N N hv− , w− i
−
N −1 X
n hv− , [w]n i
N Z X + −hv, w0 i + A(v, Bα∗ w) dt.
n=1
n=1
(3.2)
In
Since −hv, z 0 i + A(v, Bα∗ z) = hv, −z 0 + Bα∗ Azi = 0, the solution z of the dual problem (1.5) satisfies N GN (v, z) = hv− , zT i for all v ∈ C J, D(A1/2 ) . (3.3) We therefore define the DG solution Z ∈ W of (1.5) by N GN (V, Z) = hV−N , Z+ i for all V ∈ W,
(3.4)
N with Z+ = zT , and deduce the Galerkin orthogonality property
GN (V, Z − z) = 0
for all V ∈ W.
(3.5)
The following representation is the basis for our analysis of the nodal error. Theorem 3.3. If u and z are the solutions of the initial-value problem (1.1) and of the dual problem (1.5), and if U and Z are the corresponding DG solutions, then N hU− − u(tN ), zT i = GN (u − Π− u, Z − z)
for every zT ∈ H.
Proof. Taking V = U in (3.4) and v = u in (3.3) gives N N hU− − u(tN ), zT i = hU− , zT i − hu(tN ), zT i = GN (U, Z) − GN (u, z) = GN (u, Z − z),
where the last step used the Galerkin orthogonality property (2.8) of U , with X = Z. Now use the Galerkin orthogonality property (3.5) of Z, with V = Π− u. 3.3. Error in the DG solution of the dual problem. We will use the following regularity estimates. Lemma 3.4. For −1 < α < 1 and 0 < t < T , the solution z of the dual problem (1.5) satisfies kA−1 z 0 (t)k + (T − t)kA−1 z 00 (t)k ≤ C(T − t)α kzT k and (T − t)1+α kAz(t)k + kz(t)k + (T − t)kz 0 (t)k ≤ CkzT k. Proof. Define the time reversal operator Rv(t) = v(T − t). Since R∂t = −∂t R and RBα∗ = Bα R, we deduce from (1.5) that the function v = RA−1 z satisfies v 0 + Bα Av = 0
for 0 < t < T ,
with v(0) = A−1 zT .
Known results for −1 < α < 0 [19, Theorem 4.2] and 0 < α < 1 [21, Theorem 2.1 with r = 2 and µ = −2] give kv 0 (t)k + tkv 00 (t)k ≤ Ctα kAv(0)k = Ctα kzT k, implying the first estimate. Similarly [19, Theorem 4.1], the function w = Rz satisfies t1+α kAw(t)k + kw(t)k + tkw0 (t)k ≤ Ckw(0)k = CkzT k,
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
9
implying the second estimate. To investigate the DG error for the dual problem, we make the splitting A−1 (Z − z) = ζ + Θ
where ζ = A−1 (Π+ z − z) and Θ = A−1 (Z − Π+ z) ∈ W. (3.6)
Lemma 3.5. The function ζ in (3.6) satisfies α
kζkJ ≤ CtN+ k 1+α− kzT k . Proof. By (2.15) and Lemma 3.4, kζkJ is bounded by Z Z k(I − Π+ )A−1 zkJ ≤ 3 max kA−1 z 0 (t)k dt ≤ CkzT k max (tN − t)α dt. 1≤n≤N
1≤n≤N
In
In
If −1 < α < 0, then Z (1 + α) (tN − t)α dt = (tN − tn−1 )1+α − (tN − tn )1+α ≤ kn1+α , In
R whereas if 0 < α < 1, then In (tN − t)α dt ≤ kn tα N. Lemma 3.6. The function Θ in (3.6) satisfies kΘk2J
Z ≤ 8
T
0
hΘ(t), Bα∗ Aζi dt .
Proof. By (3.5), GN (V, ζ + Θ) = G(A−1 V, Z − z) = 0 for all V ∈ W, where we used the identity GN (v, A−1 w) = GN (A−1 v, w) and the fact that A−1 V ∈ W. Thus, GN (V, Θ) = −G(V, ζ)
for all V ∈ W.
n The first property in (2.12) means that ζ+ = 0 for 0 ≤ n ≤ N − 1, so the formula (3.2) shows
GN (V, ζ) =
N X
n hV−n , ζ− i+
n=1
N Z X −hV, ζ 0 i + A(V, Bα∗ ζ) dt, n=1
In
R R n n and integration by parts gives In hV, ζ 0 i dt = hV−n , ζ− i − In hV 0 , ζi dt = hV−n , ζ− i, 0 where, in the last step, we used the second property in (2.12) and the fact that V is constant on In . Thus, if we define g = −ABα∗ ζ = −Bα∗ Aζ then Z GN (V, Θ) = − 0
T
A(V, Bα∗ ζ) dt
Z
T
hV, gi dt
=
for all V ∈ W,
0
which means that Θ ∈ W is the DG solution of −θ0 + Bα∗ Aθ = g(t) for 0 < t < T , with θ(T ) = 0. The desired estimate follows by the stability of Θ, which we can prove by applying Theorem 2.1 to RΘ(t) = Θ(T − t). Lemma 3.7. If −1 < α < 1 then Z T α ∗ hV, Bα Aζi dt ≤ CtN+ k 1+α− `(tN /kN )kV kJ kzT k for all V ∈ W. 0
10
KASSEM MUSTAPHA AND WILLIAM MCLEAN
n−1 Proof. Suppose first that −1 < α < 0. Since Bα V = (B1+α V )0 and ζ+ = 0, we see using (3.1) and integrating by parts that N Z X
T
Z
hV, Bα∗ Aζi dt =
0
h(B1+α V )0 , Aζi dt =
In
n=1
N Z X
h∆n , Aζ 0 i dt,
In
n=1
where, for t ∈ In , ∆n (t) = B1+α V (tn ) − B1+α V (t) Z Z t [ω1+α (tn − s) − ω1+α (t − s)]V (s) ds + =
tn
ω1+α (tn − s)V (s) ds.
t
0
The function ω1+α is monotone decreasing whereas ω2+α is monotone increasing, so t
Z tn ω1+α (tn − s) ds k∆ (t)k ≤ kV kJ ω1+α (t − s) − ω1+α (tn − s) ds + 0 t = kV kJ ω2+α (t) − ω2+α (tn ) + 2ω2+α (tn − t) ≤ 2kV kJ ω2+α (tn − t) Z
n
R T and the Cauchy–Schwarz inequality shows that 0 hV, Bα∗ Aζi dt is bounded by 2kV kJ
NX −1
Z
ω2+α (tN − t)kAζ k dt .
Z
0
0
kAζ k dt +
ω2+α (kn ) In
n=1
IN
Since Aζ = Π+ z − z, the representation (2.14) and Lemma 3.4 imply N −1 X
Z
0
kAζ k dt ≤ C
ω2+α (kn ) In
n=1
=C
N −1 X
kn1+α
Z
n=1
In
N −1 X
Z
kn1+α
2 kz (t)k + 2 kn 0
1+
≤ Ck 1+α kzT k
Z
(tn − s)kz (s)k ds dt 0
In
In
n=1
Z
tN −1
2(tn − t) kz 0 (t)k dt kn
(T − t)−1 dt = Ck 1+α log(tN /kN )
0
and Z IN
Z 2 0 0 (T −s)kz (s)k ds dt ω2+α (tN −t)kAζ k dt ≤ C (T −t) kz (t)k+ 2 kN IN IN Z α 2kN 1+α =C (T − t)α + (T − t)kz 0 (t)k dt ≤ CkzT kkN . 2+α IN Z
0
1+α
The desired estimate follows at once. Now let 0 < α < 1. By part 2 of Lemma 3.1, Z
0
T
Z hV, Bα∗ Aζi dt ≤
T
Z kV (t)k
0
Z ≤ kV kJ
T
ωα (s − t)kAζ(s)k ds dt t
Z kAζ(s)k
0
T
s α
Z
ωα (s − t) dt ds ≤ CT kV kJ 0
T
kAζ(t)k dt. 0
11
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
The estimates (2.15) and (2.16) imply that T
Z
kAζ(t)k dt ≤ 0
N X
kn kAζkIn ≤ 4kN kzkIN + 3
n=1
N −1 X
Z kn
n=1
kz 0 (t)k dt,
In
and we know from Lemma 3.4 that kzkIN ≤ CkzT k and N −1 X
Z
Z
0
tN −1
kz (t)k dt ≤ CkkzT k
kn In
n=1
(tN − t)−1 dt = CkkzT klog(tN /kN ).
0
Hence, we arrive at the following error estimate for the dual problem. Theorem 3.8. Let z denote the solution of the dual problem (1.5), and let Z denote the DG solution defined by (3.4). Then, for −1 < α < 1, α
kA−1 (Z − z)kJ ≤ CtN+ k 1+α− `(tN /kN )kzT k. Proof. The splitting (3.6) implies that kA−1 (Z − z)kJ ≤ kζkJ + kΘkJ , and we estimate these two terms using Lemmas 3.5, 3.6 and 3.7. 4. Nodal superconvergence. With the help of Theorems 3.3 and 3.8, we are n ≈ u(tn ). Define now able to estimate the error in the approximation U− (u) = D1 + E1 + max kj Dj + 2≤j≤n
n X
2+α−
kj
Ej ,
(4.1)
j=2
where Z D1 =
kAu0 (t)k dt
Z and E1 =
t1+α− kA2 u0 (t)k dt,
(4.2)
I1
I1
with Z Dn =
kAu00 (t)k dt
Z and En =
In
kA2 u00 (t)k dt
for 2 ≤ n ≤ N .
(4.3)
In
For the rest of the paper, we assume that k1 ≤ k2 ≤ · · · ≤ kN .
(4.4)
Theorem 4.1. Let u be the solution of the initial value problem (1.1) and let U be the DG solution satisfying (2.4). Then, for 1 ≤ n ≤ N , n α+ 1+α− + kU− − u(tn )k ≤ Ctα `(tn /kn ) (u). n (1 + tn )k
(4.5)
Proof. Put η = u − Π− u and define Z Z n 0 n δ1 = hη, (z − Z) i dt and δ2 = hBα Aη, z − Zi dt. In
Since
n η−
In
= 0 for 1 ≤ n ≤ N , we see from Theorem 3.3, (3.2) and Lemma 3.1 that N hU− − u(tN ), zT i = GN (η, Z − z) =
N X n=1
δ1n + δ2n .
12
KASSEM MUSTAPHA AND WILLIAM MCLEAN
Since Z 0 is constant on In , the second property of Π− in (2.10) gives Z Z Z Z t
n 0 0 0 Aη(t), A−1 z 00 (s) ds dt, δ1 = hη, z i dt = η(t), z (t) − z (tn−1 ) dt = In
In
tn−1
In
and therefore, using Lemma 3.4, N −1 X
|δ1n |
N −1 X
≤ kAηkJ
n=1
Z
whereas N −1 X
Z kn
R =
z (t)k dt ≤ CkAηkJN kzT k
N −1 X
In
n=1
|δ1N |
kA
kn
−1 00
Z
(tN −t)α−1 dt,
kn
n=1
In
R z i dt ≤ CkAηkIN kzT k IN (tN − t)α dt. Here,
−1 0
IN
hAη, A
(tN − t)α−1 dt
In
n=1
−α−
≤ k 1+α− kN
−α−
tN −1
Z
(tN − t)α−1 dt =
0
k 1+α− kN α
α
+ 1+α− α (tα , N − kN ) ≤ CtN k
R α 1+α and likewise IN (tN − t)α dt = kN /(1 + α) ≤ CtN+ k 1+α− . By (2.15), kAηkI1 ≤ 3D1 and kAηkIn ≤ 2kn Dn for 2 ≤ n ≤ N , so N X
α |δ1n | ≤ CtN+ k 1+α− kzT k D1 + max kn Dn for −1 < α < 1. 2≤n≤N
n=1
(4.6)
Turning to δ2n , if −1 < α < 0, then [31, Lemma 2] N Z X n δ2 = n=1
tN 2
−1
hBα A η, A
0
N X −1 2+α (Z − z)i dt ≤ CkA (Z − z)kJ E1 + kn En , n=2
(4.7) but if 0 < α < 1 then |δ2n | is bounded by Z Z −1 2 −1 kA (Z − z)kIn kBα A η(t)k dt ≤ kA (Z − z)kIn In
In
Z
t
ωα (t − s)kA2 η(s)k ds,
0
so, after summing over n and reversing the order of integration, N X
−1 |δ2n | ≤ Ctα (Z − z)kJ N kA
Z
tN
kA2 η(t)k dt.
0
n=1
The integral representation (2.11) implies that Z Z Z Z t1 t1 − t skA2 u0 (s)k ds dt kA2 η(t)k ≤ kA2 u0 (s)k ds + 2 2 k1 I1 I1 I t Z s Z1 Z s 2 0 = kA u (s)k dt + 2 2(t1 − t) dt ds = 2E1 , k1 I1 I1 0 Rt PN PN and by (2.15), t1N kA2 η(t)k dt ≤ n=2 kn kA2 ηkIn ≤ n=2 2kn2 Ej . Combine this with (4.7) and applying Theorem 3.8, N X n=1
N X 2α |δ2n | ≤ CtN + k 1+α− `(tN /kN )kzT k E1 + kn2+α− En for −1 < α < 1. (4.8) n=2
13
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
Since zT ∈ H is arbitrary, the desired estimate follows from (4.6) and (4.8). To estimate the convergence rate at the nodes, we introduce some assumptions about the behaviour of the time steps, namely that, for some fixed γ ≥ 1, kn ≤ Cγ k min(1, t1−1/γ ) n
and tn ≤ Cγ tn−1
for 2 ≤ n ≤ N ,
(4.9)
with cγ k γ ≤ k1 ≤ Cγ k γ .
(4.10)
For example, these assumptions are satisfied if we put tn = (n/N )γ T
for 0 ≤ n ≤ N .
(4.11)
Recall that `(t) = max(1, | log t|). Lemma 4.2. Assume that u satisfies (1.3) and (1.4), and that the time mesh satisfies (4.9) and (4.10). Then, with γ ∗ = (2 + α− )/σ and for 1 ≤ n ≤ N , γσ 1 ≤ γ < γ∗, k , (u) ≤ CT M × k 2+α− `(tn /k1 ), γ = γ ∗ , 2+α− k , γ > γ∗. Rk Proof. The stated assumptions imply that D1 + E1 ≤ CM 0 1 tσ−1 dt ≤ CM k1σ ≤ CM k γσ , and, for 2 ≤ j ≤ n, ( Z k γσ , 1 ≤ γ < 2/σ, σ−2 2 σ−2 kj Dj ≤ CM kj t dt ≤ CM kj tj ≤ CM × σ−2/γ 2 k , γ ≥ 2/σ, tn Ij γσ 1 ≤ γ ≤ γ∗, k , 2+α − ≤ CM × k , γ ∗ ≤ γ ≤ 2/σ, σ−(2+α− )/γ 2+α− , γ ≥ 2/σ. tn k Similarly, n X j=2
2+α kj − Ej
≤M
n X
2+α kj −
j=2
≤ CM k 2+α−
Z t
σ−3−α−
dt ≤ CM k
2+α−
Z
Ij
tn
tσ−1−(2+α− )/γ dt
k1
γσ−(2+α− ) , 1 < γ < γ∗, k × log(tn /k1 ), γ = γ∗, σ−(2+α− )/γ tn , γ > γ∗.
We can now state our main result on nodal superconvergence. Theorem 4.3. Assume that the solution u of the initial value problem (1.1) satisfies (1.3) and (1.4), and that the time mesh satisfies (4.9) and (4.10). Then, with γ ∗ = (2 + α− )/σ and for 1 ≤ n ≤ N , the DG method (2.4) satisfies the nodal error bound 1+α− +γσ `(tn /kn ), 1 ≤ γ < γ ∗ , k n 3+2α − kU− − u(tn )k ≤ CT M × k `(tn /kn )2 , γ = γ∗, 3+2α− k `(tn /kn ), γ > γ∗ . Proof. The error bound follows at once from Theorem 4.1 and Lemma 4.2.
14
KASSEM MUSTAPHA AND WILLIAM MCLEAN
5. Postprocessing. We can postprocess the DG solution U to obtain a globally superconvergent solution U ] using simple Lagrange interpolation, as follows. Given a piecewise continuous function v : J → H, define Lv : J → H by linear interpolation on the first two subintervals, n−1 n Lv(t) = kn−1 [(tn − t)v− + (t − tn−1 )v− ]
for t ∈ In and n ∈ {1, 2},
(5.1)
and backward quadratic interpolation on the remaining subintervals, (t − tn−1 )(t − tn ) n−2 (t − tn−2 )(t − tn ) n−1 (t − tn−2 )(t − tn−1 ) n v − v− + v− kn−1 (kn−1 + kn ) − kn−1 kn (kn−1 + kn )kn (5.2) n for t ∈ In and n ≥ 3. Thus, (Lv)(tn ) = v− for 0 ≤ n ≤ N , and we define the postprocessed solution by Lv(t) =
U ] = LU.
(5.3)
The interpolant of the exact solution satisfies the following error bound. Lemma 5.1. Assume that there exist positive constants M and σ ] such that the regularity estimate (1.6) holds, and put γ ] = 3/σ ] . If the time mesh satisfies (4.9)– (4.10), then ( ] k γσ , 1 ≤ γ < γ], ku − LukJ ≤ CM × σ ] −3/γ 3 T k , γ ≥ γ]. Proof. If n ∈ {1, 2} and t ∈ In , then Z
t
u0 (s) ds −
(u − Lu)(t) = tn−1
t kn
Z
tn
u0 (s) ds,
tn−1
and thus, noting that t2 ≤ Cγ t1 , ku − LukIn ≤
Z tn ] ] tn 1+ ku0 (s)k ds ≤ CM tσ1 ≤ CM k γσ . kn tn−1
If n ≥ 3 and t ∈ In , then we can write the interpolation error in terms of a divided difference, (u − Lu)(t) = u[tn−2 , tn−1 , t, tn ](t − tn−2 )(t − tn−1 )(t − tn ), so ]
1 ku − LukIn ≤ 41 kn2 (kn−1 + kn ) 3! ku000 k[tn−2 ,tn ] ≤ CM kn2 (kn−1 + kn )tσn
−3
,
where, in the final step, we used (4.9). If 1 ≤ γ < 3/σ ] then, again using (4.9), ]
kn2 (kn−1 + kn )tσn
−3
]
]
]
≤ C(kt1−1/γ )γσ kn3−γσ tσn n
−3
]
]
]
= Ck γσ (kn /tn )3−γσ ≤ Ck γσ ,
but for γ ≥ 3/σ, ]
kn2 (kn−1 + kn )tσn
−3
]
≤ C(kt1−1/γ )3 tσn n
−3
]
≤ Ck 3 tσn
−3/γ
≤ CT σ
]
−3/γ 3
k .
Now consider the stability of the interpolation operator L. We see from (5.1) that 0 1 1 2 kLvkI1 ≤ max |v− |, |v− | and kLvkI2 ≤ max |v− |, |v− | .
15
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
A similar estimate holds for the subsequent subintervals provided the mesh satisfies the local quasi-uniformity condition kn ≤ Λkn−1
for 3 ≤ n ≤ N .
(5.4)
For example, our standard mesh (4.11) satisfies this condition with Λ = 2γ − 1. Lemma 5.2. If (5.4) holds, then j kLvkIn ≤ 2 + 45 Λ max |v− | for 3 ≤ n ≤ N . n−2≤j≤n
Proof. The estimate follows from (5.2) because, for t ∈ In and n ≥ 2, 1 2 1 kn |(t − tn−1 )(t − tn )| 4 kn ≤ ≤ 4 ≤ 14 Λ, kn−1 (kn−1 + kn ) kn−1 (kn−1 + kn ) kn−1 |(t − tn−2 )(t − tn )| (kn−1 + kn )kn kn ≤ =1+ ≤ 1 + Λ, kn−1 kn kn−1 kn kn−1 (kn−1 + kn )kn |(t − tn−2 )(t − tn−1 )| ≤ = 1. (kn−1 + kn )kn (kn−1 + kn )kn
Hence, the interpolant U ] is superconvergent, uniformly in t, for appropriate γ. Theorem 5.3. Suppose that the time mesh satisfies (4.9), (4.10) and (5.4), and that u satisfies (1.6). Then, with γ ] = 3/σ ] , the postprocessed solution (5.3) satisfies ( ] k γσ , 1 ≤ γ < 3/σ ] , n ] 5 kU − ukJ ≤ (2 + 4 Λ) max kU− − u(tn )k + CT M × 0≤n≤N k3 , γ ≥ 3/σ ] . Proof. Write U ] − u = L(U − u) + (Lu − u) and apply Lemmas 5.1 and 5.2. 6. Spatial discretization. 6.1. The fully discrete DG method. We denote the norm of u in H r (Ω) by kukr , and assume now that A = −∇2 in a bounded, convex or C 2 domain Ω in Rd , subject to homogeneous Dirichlet boundary conditions. Thus, if u ∈ H01 (Ω) and Au ∈ L2 (Ω), then u ∈ H 2 (Ω) and kuk2 ≤ CkAuk (kvk2 := kvkH 2 (Ω) and kvk := kvkL2 (Ω) ). Let Sh ⊆ D(A1/2 ) = H01 (Ω) denote the space of continuous, piecewiselinear functions with respect to a quasi-uniform partition of Ω into triangular or quadrilateral (or tetrahedral etc.) finite elements, with maximum diameter h. Recall that the L2 -projector Ph : L2 (Ω) → Sh and the Ritz projector Rh : H01 (Ω) → Sh are defined by hPh v, W i = hv, W i and A(Rh v, W ) = A(v, W )
for all W ∈ Sh ,
(6.1)
and that the latter has the quasi-optimal approximation property kv − Rh vk + hk∇(v − Rh v)k ≤ Ch2 kvk2
for v ∈ H01 (Ω) ∩ H 2 (Ω).
(6.2)
Let W(Sh ) denote the space of piecewise linear functions U : J → Sh (so U is continuous in space, but may be discontinuous in time). We define the fully discrete DG solution Uh ∈ W(Sh ) by requiring (2.4) to hold for every X ∈ W(Sh ). Equivalently, cf. (2.6), Z tN 0 0 GN (Uh , X) = hUh− , X+ i+ hf (t), X(t)i dt for all X ∈ W(Sh ), (6.3) 0
16
KASSEM MUSTAPHA AND WILLIAM MCLEAN
0 where, for simplicity, we choose Uh− = (Uh )0− = Ph u0 . In view of (2.7), the Galerkin orthogonality property (2.8) now takes the form 0 GN (Uh − u, X) = hPh u0 − u0 , X+ i=0
for all X ∈ W(Sh ).
(6.4)
Similarly, the fully discrete DG solution Zh ∈ W(Sh ) for the dual problem (1.5) is defined by N GN (V, Zh ) = hV−N , Zh+ i for all V ∈ W(Sh ),
N with Zh+ = Ph zT ,
(6.5)
and, since z satisfies (3.3), GN (V, Zh − z) = hV−N , Ph zT − zT i = 0
for all V ∈ W(Sh ).
(6.6)
Theorem 3.3 generalizes as follows. Theorem 6.1. If u and z are the solutions of the initial value problem (1.1) and the dual problem (1.5), and if Uh and Zh are the corresponding fully discrete DG solutions satisfying (6.3) and (6.5), then
N Uh− − u(tN ), zT = GN (u − Π− Rh u, Zh − z) for every zT ∈ H. Proof. By taking V = Uh in (6.5) we find that N N N N hUh− , zT i = hUh− , Ph zT i = hUh− , Zh+ i = GN (Uh , Zh ),
and taking v = u in (3.3) we have hu(tN ), zT i = GN (u, z), so
N Uh− − u(tN ), zT = GN (Uh , Zh ) − GN (u, z) = GN (Uh − u, Zh ) + GN (u, Zh − z) = GN (u, Zh − z), where the final step used (6.4) with X = Zh . Since GN (u, Zh − z) = GN (u − Π− Rh u, Zh − z) + GN (Π− Rh u, Zh − z), the result follows after putting V = Π− Rh u in (6.6). 6.2. Error in the fully discrete DG solution of the dual problem. We modify the splitting (3.6), by writing A−1 (Zh − z) = ζ + Ψ + Φ where ζ = A−1 (Π+ z − z),
Ψ = A−1 Π+ (Ph z − z),
Φ = A−1 (Zh − Π+ Ph z) ∈ W(Sh ).
Theorem 3.8 generalizes as follows. Theorem 6.2. Let z denote the solution of the dual problem (1.5), and let Zh denote the fully discrete DG solution defined by (6.5). Then, for −1 < α < 1, α kA−1 (Zh − z)kJ ≤ C tN+ k 1+α− + h2 `(tN /kN )kzT k. Proof. We already estimated kζk in Lemma 3.5. To estimate Ψ, observe that since A−1 commutes with Π+ and since ARh = Ph A (implying A−1 Ph = Rh A−1 ), Ψ = Π+ (Rh − I)A−1 z.
(6.7)
17
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
Using (2.16), the error bound (6.2) for the Ritz projection, H 2 -regularity for A and Lemma 3.4, we find that kΨkIn ≤ 3k(Rh − I)A−1 zkIn ≤ Ch2 kA−1 z(t)k2 ≤ Ch2 kz(t)k ≤ Ch2 kzT k.
(6.8)
To estimate Φ, let V ∈ W(Sh ) and observe that since A−1 V = A−1 Ph V = Rh A−1 V , GN (V, ζ + Ψ + Φ) = GN (A−1 V, Zh − z) = GN (Rh A−1 V, Zh − z) = 0,
(6.9)
where we used (6.6) with V replaced by Rh A−1 V . From the proof of Lemma 3.6, Z T hV, Bα∗ Aζi dt, (6.10) GN (V, ζ) = − 0
and by (3.2), GN (V, Ψ) =
hV−N , ΨN −i
−
N −1 X
hV−n , [Ψ]n i
N Z X + −hV, Ψ0 i + A(V, Bα∗ Ψ) dt. (6.11)
n=1
n=1
In
Since A(V, Bα∗ Ψ) = A V, Bα∗ A−1 Π+ (Ph − I)z) = hV, (Ph − I)Bα∗ Π+ zi = 0, |GN (V, Ψ)| ≤ kV kJ
kΨN −k
+
N −1 X
n
k[Ψ] k +
n=1
N Z X n=1
kΨ k dt . 0
In
2 By (6.8), kΨN − k ≤ kΨkIN ≤ Ch kzT k. Using (6.7), (2.16) and (6.2) we have Z Z n −1 0 2 k[Ψ] k ≤ k(Rh − I)A z (t)k dt ≤ Ch kA−1 z 0 (t)k2 dt, In
In
so H 2 -regularity for A and the estimate (T − t)kz 0 (t)k ≤ CkzT k from Lemma 3.4 give N −1 X
k[Ψ]n k ≤ Ch2 kzT k
Z
tN −1
(T − t)−1 dt = Ch2 kzT k log(tN /kN ).
0
n=1
Using (6.7), (2.16), (6.2) and Lemma 3.4, we find that Z Z Z 2 Ch2 tn − t kΨ0 (t)k dt ≤ (tn − t)k(Rh − I)A−1 z 0 (t)k dt ≤ kzT k dt, kn In kn In In T − t so N Z X n=1
In
Z kΨ k dt ≤ Ch kzT k 0
2
0
tN −1
−1
(T − t)
dt +
−1 kN
Z
dt ,
IN
and we conclude that |GN (V, Ψ)| ≤ Ch2 `(tN /kN )kV kJ kzT k. Therefore, by (6.9), (6.10) and Lemma 3.7, α |GN (V, Φ)| = |GN (V, ζ)+GN (V, Ψ)| ≤ C tN+ k 1+α− +h2 `(tN /kN )kV kJ kzT k. (6.12) Fix n with 1 ≤ n ≤ N , and define V ∈ W(Sh ) by ( 0, if t ∈ Ij for 1 ≤ j ≤ n − 1, V (t) = Φ(t), if t ∈ Ij for n ≤ j ≤ N .
18
KASSEM MUSTAPHA AND WILLIAM MCLEAN
In view of (2.9), GN (V, Φ) ≥ mate (6.12) gives
1 N 2 2 kΦ− k
2 + 12 kΦn−1 + k +
1 2
PN −1 j=n
k[Φ]j k2 , so the esti-
α+ 1+α− 2 n 2 + h2 `(tN /kN )kΦk(tn−1 ,T ) kzT k kΦn−1 + k + k[Φ] k ≤ C tN k for 1 ≤ n ≤ N − 1, whereas α+ 1+α− −1 2 2 + h2 `(tN /kN )kΦk(tN −1 ,T ) kzT k. kΦN k + kΦN − k ≤ C tN k + n Furthermore, kΦkIn = max kΦn−1 + k, kΦ− k because Φ is piecewise linear in t, and 2 n 2 n 2 kΦn− k ≤ kΦn+ k + k[Φ]n k, implying that kΦk2In ≤ kΦn−1 + k + kΦ+ k + k[Φ] k . By letting n∗ = argmax1≤n≤N kΦkIn , we see that α kΦk2J = kΦk2In∗ ≤ C tN+ k 1+α− + h2 `(tN /kN )kΦkJ kzT k, giving the desired bound for kΦkJ . 6.3. Fully-discrete error . As claimed in the Introduction, we have the following error bound for Uh ; recall that γ ∗ = (2 + α− )/σ whereas γ ] = 3/σ ] . Theorem 6.3. Assume that the solution u of the initial value problem (1.1) satisfies (1.3) and (1.4), and that the time mesh satisfies assumptions (4.9) and (4.10). Then for 0 ≤ n ≤ N the fully discrete DG solution Uh ∈ W(Sh ) satisfies 1+α− +γσ `(tn /kn ), 1 ≤ γ < γ ∗ , k n 2 kUh− − u(tn )k ≤ CT M h + CT M × k 3+2α− `(tn /kn )2 , γ = γ∗, 3+2α− k `(tn /kn ), γ > γ∗. Moreover, if the additional regularity property (1.6) holds, then the fully discrete, postprocessed solution Uh] = LUh satisfies ( ] k γσ , 1 ≤ γ < γ ] , ] n 5 kUh − ukJ ≤ (2 + 4 Λ) max kUh− − u(tn )k + CT M × 0≤n≤N k3 , γ ≥ γ]. Proof. In view of Lemma 4.2, it suffices to show (cf. Theorem 4.1) that n kUh− − u(tn )k ≤ CT k 1+α− + h2 `(tn /kn )(u) + CT M h2 . Put ξ = u − Rh u and η = u − Π− u so that u − Π− Rh u = η + Π− ξ and thus, N by Theorem 6.1, hUh− − u(tN ), zT i = GN (η, Zh − z) + GN (Π− ξ, Zh − z). Using Theorem 6.2 in place of Theorem 3.8, we can show as in the proof of Theorem 4.1 that |GN (η, Zh − z)| ≤ CT kzT k k 1+α− + h2 `(tn /kn ) (u). By (2.5), GN (Π− ξ, Zh ) equals
−
(Π
0 ξ)0+ , Zh+
+
N −1 X
N
− n n X [Π ξ] , Zh+ +
n=1
n=1
Z
(Π− ξ)0 , Zh + A(Bα Π− ξ, Zh ) dt,
In
and, since Bα commutes with the Ritz projector Rh , the definition (6.1) of Rh implies that A(Bα Π− ξ, Zh ) = A(Bα Π− u, Zh ) − A(Rh Bα Π− u, Zh ) = 0. Integrating by parts, n and noting that (Π− ξ)n− = ξ− = ξ(tn ), Z Z
− 0
n n − n−1 n−1
− (Π ξ) , Zh dt = ξ− , Zh− − (Π ξ)+ , Zh+ − Π ξ, Zh0 dt, In
In
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
19
and, because Zh0 is constant on In , the orthogonality property of Π− in (2.10) gives Z Z Z
− n−1 n−1 0 0 n n Π ξ, Zh dt = hξ, Zh i dt = hξ− , Zh− i − hξ+ , Zh+ i − hξ 0 , Zh i dt. In
In
In
Thus, Z In
n−1 n−1 (Π− ξ)0 , Zh dt = ξ+ − (Π− ξ)n−1 + + , Zh+
Z
hξ 0 , Zh i dt
In
and so −1 N
NX
n n X 0 GN (Π− ξ, Zh ) = ξ(0), Zh+ + [ξ] , Zh+ + n=1
Z
hξ 0 , Zh i dt.
In
n=1
Using (3.3) with v = Π− ξ, and noting that [ξ]n = 0, we obtain N
X 0 GN (Π ξ, Zh − z) = ξ(0), Zh+ − ξ(T ), zT +
Z
−
hξ 0 , Zh i dt.
In
n=1
Stability of the fully discrete dual problem, kZh kJ ≤ CkzT k, follows from (2.9), giving Z T 0 GN (Π− ξ, Zh − z) ≤ CkzT k kξ(0)k + kξ(tN )k + kξ k dt 0
Z 2 ≤ Ch kzT k ku0 k2 + ku(T )k2 +
T
ku (t)k2 dt , 0
0 n where we used the error bound (6.2) for the Ritz projector. The error bound for Uh− ] 2 follows using H -regularity of A and the assumption (1.3), and the bound for Uh is proved in the same way as for U ] in Theorem 5.3.
7. Numerical results. We present a series of numerical tests using a model problem in one space dimension, of the form (1.1) with Au = −uxx , Ω = (0, 1), [0, T ] = [0, 1], and homogeneous Dirichlet (absorbing) boundary conditions. These tests reveal faster than expected convergence when α < 0, and that our regularity assumptions are more restrictive than is needed in practice. We apply the fully discrete DG method defined in Section 6.1, employing a time mesh of the form (4.11), for various choices of the mesh grading parameter γ ≥ 1, and a uniform spatial mesh consisting of M subintervals, each of length h = 1/M . We choose M = dN 3/2 e unless indicated otherwise, and always ensure that the spatial error is negligible compared to the error from the time discretization. √ 7.1. Nodal errors. We choose u0 (x) = φ1 (x)/ 2 = sin(πx) for the initial data and f (t, x) = (α + 1) tα sin(πx) for the inhomogeneous term, so that Γ(α + 2) Γ(α + 2) 2 α+1 u(t) = Eα+1 (−π t ) 1− + sin(πx), (7.1) π2 π2 P∞ where the Mittag–Leffler function is given by Eν (t) = p=0 tp /Γ(1 + νp). One easily verifies that the regularity conditions (1.3), (1.4) and (1.6) hold for σ = 1 + α = σ] ,
(7.2)
20
KASSEM MUSTAPHA AND WILLIAM MCLEAN Table 7.1 The left nodal error for Example 1 when α = −0.6, and so γ ∗ = 7/2 = 3.5.
N 20 40 80 160 320 640 ρ
N =M γ=1 6.81e-03 6.02e-03 0.177 5.10e-03 0.238 4.07e-03 0.325 3.06e-03 0.415 2.16e-03 0.498 0.800
γ=3 8.68e-04 2.41e-04 1.848 5.99e-05 2.008 1.38e-05 2.123
M = dN 3/2 e γ = 7/2 5.98e-04 1.50e-04 2.066 3.42e-05 2.171 7.37e-06 2.236
γ=4 4.76e-04 1.14e-04 2.141 2.51e-05 2.216 5.30e-06 2.268
1.600
1.800
1.800
Table 7.2 Left nodal errors when α = −0.25, and so γ ∗ = 7/3.
N 20 40 80 160 ρ
γ=1 5.37e-03 1.83e-03 1.556 5.58e-04 1.712 1.66e-04 1.752 1.500
γ=2 2.57e-04 4.27e-05 2.588 6.78e-06 2.653 1.05e-06 2.685 2.250
γ = 7/3 2.18e-04 3.55e-05 2.712 5.58e-06 2.717 8.62e-07 2.719 2.500
γ=3 2.05e-04 3.34e-05 2.714 5.27e-06 2.711 8.13e-07 2.721 2.500
Table 7.3 Right nodal errors when α = −0.25.
N 20 40 80 160
γ=1 1.31e-01 7.74e-02 0.809 4.42e-02 0.809 2.50e-02 0.821
γ = 7/3 6.69e-03 1.73e-03 2.017 4.91e-04 1.853 1.44e-04 1.789
γ=3 4.94e-03 1.26e-03 2.039 3.20e-04 2.016 8.06e-05 2.008
γ=4 4.82e-03 1.23e-03 2.035 3.14e-04 2.013 7.91e-05 2.006
because Au(t) = π 2 u(t) and A2 u(t) = π 4 u(t). n = limt→t± Uh (t), we use the terminology With Uh± n n left nodal error = max ||Uh− − u(tn )||, 1≤n≤N
n right nodal error = max ||Uh+ − u(tn )||, 1≤n≤N
and expect the former, but not the latter, to exhibit superconvergence. By Theorem 6.3, the left nodal error is of order k ρ (ignoring the logarithmic factors) where ( 1 + α− + γσ, 1 ≤ γ < γ ∗ = (2 + α− )/σ, ρ= (7.3) 3 + 2α− , γ ≥ γ∗, provided h2 k ρ . Tables 7.1 and 7.2 show the left nodal errors when α = −0.6 and −0.25, respectively, for various choices of N and γ. We observe that the convergence rates are faster than predicted by our theory, for σ given by (7.2), and appear to be O(k 3+α ) when γ ≥ γ ∗ . Table 7.3 shows the right nodal errors when α = −0.25, which never exhibit better than O(k 2 ) convergence, and thus are not superconvergent. Tables 7.4 and 7.5 show the left and right nodal errors, respectively, when α = +0.25, and Table 7.6 shows the left nodal errors when α = +0.75. For these positive
21
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD Table 7.4 Left nodal errors when α = +0.25, and so γ ∗ = 8/5 = 1.6.
N 20 40 80 160 ρ
γ=1 4.98e-05 1.22e-05 2.103 2.57e-06 2.288 4.91e-07 2.408 2.250
γ = 1.3 4.33e-05 5.20e-06 3.170 6.23e-07 3.116 7.54e-08 3.074 2.625
γ = 8/5 4.45e-05 5.38e-06 3.157 6.52e-07 3.102 7.94e-08 3.065 3.000
γ = 1.9 5.11e-05 6.22e-06 3.146 7.60e-07 3.087 9.34e-08 3.053 3.000
Table 7.5 Right nodal errors when α = +0.25.
N 20 40 80 160
γ=1 1.00e-02 5.22e-03 0.973 2.39e-03 1.149 1.04e-03 1.211
γ = 1.3 4.17e-03 1.46e-03 1.572 4.84e-04 1.619 1.58e-04 1.628
γ = 8/5 3.04e-03 7.67e-04 2.059 1.93e-04 2.027 4.84e-05 2.015
γ = 1.9 3.27e-03 8.33e-04 2.044 2.09e-04 2.027 5.25e-05 2.015
Table 7.6 Left nodal errors when α = +0.75, and so γ ∗ = 8/7 ≈ 1.143.
N 20 40 80 160 ρ
γ=1 5.85e-05 7.51e-06 3.067 9.46e-07 3.042 1.19e-07 3.018 2.750
γ = 1.1 5.48e-05 6.99e-06 3.080 8.77e-07 3.047 1.10e-07 2.995 2.925
γ = 8/7 5.42e-05 6.90e-06 3.081 8.65e-07 3.048
γ = 1.3 5.34e-05 7.00e-06 3.036 8.95e-07 3.022
3.000
3.000
values of α, the left nodal errors appear to be O(k 3 ) when γ is sufficiently large, whereas the right nodal errors are again at best O(k 2 ). 7.2. Uniform error after post-processing. We introduce a finer mesh G N,m = { tj−1 + `kj /m : j = 1, 2, . . . , N and ` = 0, 1, . . . , m },
(7.4)
and define the discrete maximum norm kvkJ,m = maxt∈G N,m kv(t)k, so that, for sufficiently large values of m, kUh] − ukJ,m approximates the uniform error kUh] − ukJ , ] which, by Theorem 6.3, is of order k ρ where ( min(ρ, γσ ] ), 1 ≤ γ < γ ] = 3/σ ] , ] ρ = min(ρ, 3), γ ≥ γ]. Tables 7.7 and 7.8 show kUh] − ukJ,12 when α = −0.25 and α = +0.25, respectively. Comparing these results with those of Tables 7.2–7.5 confirms that the full rate of nodal superconvergence is reflected in the uniform error after postprocessing, but only with a stronger mesh grading, since γ ] > γ ∗ . 7.3. A less regular solution. Now let u0 (x) = x(1 − x) and f ≡ 0. Thus, Au0 = 2 does not vanish at 0 and 1, so u0 ∈ / D(A2 ) [19, Equation (3.2)]. Separation of variables yields a series representation u(x, t) = 8
∞ X n=0
ωn−3 sin(ωn x)E1+α (−ωn2 t1+α )
with ωn = (2n + 1)π,
(7.5)
22
KASSEM MUSTAPHA AND WILLIAM MCLEAN Table 7.7 Uniform errors after postprocessing when α = −0.25, so γ ∗ = 7/3 and γ ] = 4.
N 20 40 80 160 320 ρ]
γ=2 8.81e-03 2.95e-03 1.635 1.02e-03 1.559 3.58e-04 1.525 1.26e-04 1.505 1.500
γ = 7/3 4.98e-03 1.30e-03 2.006 3.67e-04 1.859 1.07e-04 1.791 3.17e-05 1.755 1.750
γ=3 2.76e-03 5.48e-04 2.418 1.14e-04 2.308 2.39e-05 2.274 5.01e-06 2.253 2.250
γ=4 1.75e-03 2.14e-04 3.140 2.71e-05 3.032 3.51e-06 2.976 4.60e-07 2.931 2.500
Table 7.8 Uniform errors after postprocessing when α = +0.25, and so γ ∗ = 8/5 = 1.6 and γ ] = 12/5 = 2.4.
N 20 40 80 160 ρ]
γ = 1.3 3.11e-03 1.07e-03 1.590 3.55e-04 1.625 1.16e-04 1.630 1.625
γ = 8/5 1.08e-03 3.17e-04 1.827 8.28e-05 1.972 2.09e-05 2.002 2.000
γ = 1.9 6.49e-04 1.35e-04 2.353 2.63e-05 2.398 5.08e-06 2.393 2.375
γ = 12/5 6.43e-04 7.89e-05 3.135 9.63e-06 3.089 1.18e-06 3.052 3.000
and it is not difficult to verify directly that (1.3) holds. In fact, differentiating (7.5), ∂tj Au(x, t) = ∂tj uxx = −8
∞ X
ωn−1 sin(ωn x)
n=0
dj E1+α (−ωn2 t1+α ), dtj
and thus by Parseval’s identity, k∂tj Au(t)k2 = 32
∞ X
ωn−2
n=0
2 dj 2 1+α E (−ω t ) . 1+α n dtj
P∞ Since |E1+α (−ωn2 t1+α )| ≤ C for all t > 0 it follows that kAu(t)k2 ≤ C n=0 ωn−2 < ∞, and the Mittag–Leffler function satisfies [19, Theorem 4.2] j d −(1+α)µ−j −2µ 2 1+α ωn for j ∈ {1, 2, 3, . . .} and |µ| ≤ 1, (7.6) E (−ω t ) n dtj 1+α ≤ Ct so j
t
2 k∂tj Au(t)k
−2µ(1+α)
≤ Ct
∞ X
ωn−4µ−2
n=0
C t−µ(1+α) ≤ 1 + 4µ
2 for − 14 < µ ≤ 1.
Thus, the regularity conditions (1.3) hold for σ < 41 (1+α) and M = C(1+α−4σ)−1/2 . Similarly, ∂tj A2 u(x, t) = ∂tj uxxxx = 8
∞ X
ωn sin(ωn x)
n=0
and
k∂tj A2 u(t)k2
= 32
P∞
2 n=0 ωn
dj E1+α (−ωn2 t1+α ) dtj
2
dj 2 1+α ) dtj E1+α (−ωn t
, so
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
23
Table 7.9 Left nodal errors for the less regular solution (7.5).
N 20 40 80 160 240 ρ
γ=1 1.84e-03 1.68e-03 0.129 1.44e-03 0.223 1.15e-03 0.320 9.84e-04 0.393 0.500
N 20 80 160 ρ
γ=1 2.01e-03 3.90e-04 1.143 2.21e-04 0.821 0.875
N 20 40 80 160 ρ
γ=1 2.10e-04 6.77e-05 1.632 2.19e-05 1.636 7.11e-06 1.625 1.025
N 20 40 80 160
γ=1 7.37e-05 1.83e-05 2.011 4.56e-06 2.003 1.14e-06 2.001
j
t
2 k∂tj A2 u(t)k
α = −0.6, γ ∗ = 14 γ=3 γ=5 2.47e-04 1.09e-04 6.99e-05 1.818 2.54e-05 2.103 2.26e-05 1.628 5.52e-06 2.202 8.06e-06 1.489 1.16e-06 2.250 4.38e-06 1.500 4.60e-07 2.282 0.700 0.900 ∗ α = −0.3, γ = 68/7 ≈ 9.71 γ=2 γ=3 1.08e-04 6.39e-05 9.33e-06 1.758 1.81e-06 2.596 2.77e-06 1.753 2.92e-07 2.632 1.050 1.225 α = +0.3, γ ∗ = 80 γ = 1.5 γ = 1.75 2.08e-05 1.21e-05 3.61e-06 2.527 1.61e-06 2.904 6.43e-07 2.486 2.13e-07 2.917 1.17e-07 2.461 2.80e-08 2.930 1.038 1.044 α = +0.6 γ = 1.25 γ = 1.5 1.88e-05 1.30e-05 3.12e-06 2.590 1.87e-06 2.801 5.31e-07 2.554 2.60e-07 2.845 9.21e-08 2.528 3.64e-08 2.840
−2(1+α)µ
≤ Ct
∞ X n=0
ωn2−4µ ≤
C(t−(1+α)µ )2 4µ − 3
γ=6 1.05e-04 2.44e-05 2.111 5.29e-06 2.203 1.11e-06 2.247 2.30e-07 2.278 1.000 γ = 3.25 6.39e-05 1.82e-06 2.595 2.94e-07 2.632 1.269 γ=2 1.23e-05 1.57e-06 2.966 1.99e-07 2.983 2.53e-08 2.972 1.050 γ=2 2.19e-05 3.14e-06 2.803 4.22e-07 2.894 5.51e-08 2.938
for
3 4
< µ ≤ 1.
Thus, (1.4) holds for σ < 1 + α− − 34 (1 + α), but σ > 0 only if −1 < α < 31 . Table 7.9 shows the left nodal errors for α = −0.6, −0.3, +0.3 and +0.6. To compute the error, we evaluated u by truncating the Fourier series (7.5) after 160 terms. In the first three cases, α < 13 and we show the values of ρ given by (7.3) with σ = 1 + α− − 43 (1 + α). Evidently, the theoretical convergence rates are even more pessimistic than in our previous example (7.1), owing to the smaller value of σ, and the values of γ ∗ become unrealistically large. However, in practice, we observe that the errors are superconvergent of order k 3+α− for moderate γ, even in the case α = +0.6 when the regularity assumption (1.4) fails completely. 8. Concluding remarks. We have analysed a piecewise-linear DG method for the time discretization of (1.1) — a fractional diffusion (−1 < α < 0) or wave (0 < α < 1) equation — and proved superconvergence at the nodes, generalizing a known result for the classical heat equation. Numerical experiments indicate that our theoretical error bounds are sharp if α > 0, but not if α < 0. For generic regular
24
KASSEM MUSTAPHA AND WILLIAM MCLEAN
data u0 and f , derivatives of the exact solution are singular as t → 0, but nevertheless by employing non-uniform time steps we achieve a high convergence rate of O(k 3+α− ). After postprocessing the solution, the same high accuracy is achieved for all t, not just at the nodes. We have also proved that the additional error arising from a spatial discretization by continuous piecewise-linear finite elements is O(h2 ). REFERENCES [1] C.-M. Chen, F. Liu, V. Anh, and I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput., 32 (2010), pp. 1740–1760. , Numerical methods for solving a two-dimensional variable-order anomalous subdiffu[2] sion equation, Math. Comp., 81 (2012), pp. 345–366. [3] E. Cuesta, C. Lubich, and C. Palencia, Convolution quadrature time discretization of fractional diffusion-wave equations, Math. Comp., 75 (2006), pp. 673–696. [4] E. Cuesta and C. Palencia, A fractional trapezoidal rule for integro-differential equations of fractional order in Banach spaces, Appl. Numer. Math., 45 (2003), pp. 139–159. , A numerical method for an integro-differential equation with memory in Banach spaces: [5] qualitative properties, SIAM J. Numer. Anal., 41 (2003), pp. 1232–1241. [6] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys., 228 (2009), pp. 7792–7804. [7] M. Cui, Compact alternating direction implicit method for two-dimensional time fractional diffusion equation, J. Comput. Phys., 231 (2012), pp. 2621–2633. [8] , Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algor., (Published online: 2012). ´e, Time discretization of parabolic problems by the [9] K. Eriksson, C. Johnson, and V. Thome discontinuous Galerkin method, M2AN Math. Model. Numer. Anal., 19 (1985), pp. 611– 643. [10] A. Hanyga, Wave propagation in media with singular memory, Math. Comput. Modelling., 34 (2001), pp. 1399–1421. [11] B. Jin, R. Lazarov, and Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations. Preprint, arXiv:1204.38884v1. [12] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and Applications of Fractional Differential Equations, vol. 204 of North-Holland Mathematics Studies, North–Holland, Amsterdam, 2006. [13] T. A. M. Langlands and B. I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys., 205 (2005), pp. 719–736. [14] F. Liu, C. Yang, and K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math., 231 (2009), pp. 160–176. ´ pez-Ferna ´ ndez and C. Palencia, On the numerical inversion of the Laplace transform [15] M. Lo of certain holomorphic functions, Appl. Numer. Math., 51 (2004), pp. 289–303. ´ pez-Fernandez, C. Palencia, and A. Scha ¨ dle, A spectral order method for inverting [16] M. Lo sectorial Laplace transforms, SIAM J. Numer. Anal., 44 (2006), pp. 1332–1350. [17] F. Mainardi and P. Paradisi, Fractional diffusive waves, J. Comput. Acoustics, 9 (2001), pp. 1417–1436. [18] A. M. Mathai, R. K. Saxena, and H. J. Haubold, The H-Function: Theory and Applications, Springer, 2010. [19] W. McLean, Regularity of solutions to a time-fractional diffusion equation, ANZIAM J., 52 (2010), pp. 123–138. [20] , Fast summation by interval clustering for an evolution equation with memory. Preprint, arXiv:1203.4032v1, 2012. [21] W. McLean and K. Mustapha, A second-order accurate numerical method for a fractional wave equation, Numer. Math., 105 (2007), pp. 481–510. [22] , Convergence analysis of a discontinuous Galerkin method for a sub-diffusion equation, Numer. Algor., 52 (2009), pp. 69–88. ´e, Numerical solution of an evolution equation with a positive-type [23] W. McLean and V. Thome memory term, ANZIAM J., 35 (1993), pp. 23–70. ´e, Numerical solution via Laplace transforms of a fractional order [24] W. McLean and V. Thome
SUPERCONVERGENCE OF A DISCONTINUOUS GALERKIN METHOD
25
evolution equation, J. Integral Equations Appl., 22 (2010), pp. 57–94. ´e, and L. B. Wahlbin, Discretization with variable time steps of an [25] W. McLean, V. Thome evolution equation with a positive-type memory term, J. Comput. Appl. Math., 69 (1996), pp. 49–69. [26] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Reports, 339 (2000), pp. 1–77. , The restaurant at the end of the random walk: Recent developments in the description [27] of anomalous transport by fractional dynamics, J. Phys. A, 37 (2004), pp. R161–R208. [28] K. Mustapha, An implicit finite difference time-stepping method for a sub-diffusion equation, with spatial discretization by finite elements, IMA Journal Numer. Anal., 31 (2011), pp. 719–739. [29] K. Mustapha and J. AlMuttawa, A finite difference method for an anomalous sub-diffusion equation, theory and applications, Numer. Algor., (2012). [30] K. Mustapha and W. McLean, Discontinuous Galerkin method for an evolution equation with a memory term of positive type, Math. Comp., 78 (2009), pp. 1975–1995. , Piecewise-linear, discontinuous Galerkin method for a fractional diffusion equation, [31] Numer. Algor., 56 (2011), pp. 159–184. , Uniform convergence for a discontinuous Galerkin, time stepping method applied to a [32] fractional diffusion equation, IMA J. Numer. Anal., 32 (2012), pp. 906–925. [33] Y. nan Zhang and Z. zhong Sun, Alternating direction implicit schemes for the twodimensional fractional sub-diffusion equation, J. Comput. Phys., 230 (2011), pp. 8713– 8728. [34] I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, 1999. ¨ ss, Evolutionary Integral Equations and Applications, vol. 87 of Monographs in Mathe[35] J. Pru matics, Birkh¨ auser, Basel, 1993. [36] J. M. Sanz-Serna, A numerical method for a partial integro-differential equation, SIAM J. Numer. Anal., 25 (1988), pp. 319–327. ¨ dle, M. Lo ´ pez-Ferna ´ ndez, and C. Lubich, Fast and oblivious convolution quadra[37] A. Scha ture, SIAM J. Sci. Comput., 28 (2006), pp. 421–438. ´ ndez, and R. J. Cherry, Anoma[38] P. R. Smith, I. E. G. Morrison, K. M. Wilson, N. Ferna lous diffusion of major histocompatability complex class I molecules on hela cells determined by single particle tracking, Biophys. J., 76 (1999), pp. 3331–3344. [39] I. Sokolov and J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos, 15 (2005), p. 026103. [40] V. E. Tarasov, Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media (Nonlinear Physical Science), Springer, 2010. [41] H. Wang and K. Wang, An O(N log2 N ) alternating-direction finite difference method for two-dimensional fractional diffusion equations, J. Comput. Phys., 230 (2011), pp. 7830– 7839. [42] S. B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys., 216 (2006), pp. 264–274. [43] S. B. Yuste and L. Acedo, An explicit finite difference method and a new von Neumanntype stability analysis for fractional diffusion equations, SIAM J. Numer. Anal., 42 (2005), pp. 1862–1874. [44] P. Zhuang, F. Liu, V. Anh, and I. Turner, New solution and analytical techniques of the implicit numerical methods for the anomalous sub-diffusion equation, SIAM J. Numer. Anal., 46 (2008), pp. 1079–1095. [45] , Stability and convergence of an implicit numerical method for the nonlinear fractional reaction-subdiffusion process, IMA J. Appl. Math., 74 (2009), pp. 645–667.
