Computers and Mathematics with Applications

Contents lists available atScienceDirect

Computers and Mathematics with Applications

journal homepage:www.elsevier.com/locate/camwa

Finite element methods for fractional-order diffusion problems with optimal convergence order

Gábor Maros

^b

, Ferenc Izsák

^a,∗

aDepartment of Applied Analysis and Computational Mathematics & MTA ELTE NumNet Research Group, Eötvös Loránd University, Pázmány P. stny. 1C, 1117 Budapest, Hungary

bDepartment of Applied Analysis and Computational Mathematics, Eötvös Loránd University, Pázmány P. stny.

1C, 1117 Budapest, Hungary

a r t i c l e i n f o

Article history:

Received 27 February 2020

Received in revised form 25 June 2020 Accepted 14 September 2020 Available online 25 September 2020

Keywords:

Fractional order diffusion Matrix transformation method Finite element method Error estimation

a b s t r a c t

A convergence result is stated for the numerical solution of space-fractional diffusion problems. For the spatial discretization, an arbitrary family of finite elements can be used combined with the matrix transformation technique. The analysis covers the application of the implicit Euler method for time integration to ensure unconditional stability. The spatial convergence rate does not depend on the fractional power of the Laplacian operator. An efficient numerical implementation is developed avoiding the direct computation of matrix powers.

©2020 Elsevier Ltd. All rights reserved.

1. Introduction

Many observations confirmed the presence of the fractional diffusion in the natural sciences. Tracking the motion of food-seeking animals, the presence of ‘‘anomalous’’ diffusion was reported earlier [1,2]. In concrete terms, assuming diffusion, the flightx(t) of the individuals during timetobeys a normal distribution. Consequently, the flight-length mean

⟨|

x(t)

|⟩

of the individuals should be proportional with

√

t.

Due to accurate measurements in the last decade, the dynamics of individual molecules could also be observed: in many cases, it exhibits a subdiffusive behavior such that instead of the above proportionality, the linear relation

⟨|

x(t)

|⟩ ∼ √

t^α was detected. See a detailed overview of these measurements in [3]. A possible explanation of this dynamics can be found in [4].

A number of different mathematical models have been suggested to simulate this dynamics. Among the PDE models, in case of homogeneous Dirichlet boundary conditions, a possible choice is taking the fractional Laplacian on the entire space R^d(see [5]) and restricting to functions which are identically zero outside of the computational domain. The numerical analysis and implementation for the corresponding elliptic problems can be found in [6] and [7], respectively. In the one-dimensional case, an interesting, physically motivated approach is analyzed in [8].

Another conventional choice, which is used in the present work, is the spatial differential operator

−

(

−

_∆D)^α, where

∆^Ddenotes the Dirichlet Laplacian.

For a systematic comparison of the different approaches, we refer to [9,10] and [11].

For the numerical solution of problems with

−

(

−

_∆D)^α, the fractional power can be applied also at the discrete level, approximating (

−

_∆D)^αwith the power of the discretization of

−

_∆D. This was first observed in [12] and analyzed in [13]

∗ Corresponding author.

E-mail addresses: magabor@caesar.elte.hu(G. Maros),izsakf@caesar.elte.hu(F. Izsák).

https://doi.org/10.1016/j.camwa.2020.09.006 0898-1221/©2020 Elsevier Ltd. All rights reserved.

for the corresponding elliptic problems. A similar analysis – based on a general integral form of operator powers – was performed in [14] using an efficient numerical integration technique.

Regarding the full discretization with an implicit time stepping and efficient implementation issues, three main approaches were proposed in the literature. In [15], the authors propose a direct approximation of the solution operator (which is a matrix-power exponential) based on the elliptic theory in [14]. The corresponding contour integral is recasted as a real improper integral and approximated using an exponentially convergent quadrature.

In the works [16] and [17], explicit and implicit time discretizations are used. In [17], again using an integral representation, an efficient simple quadrature is proposed to approximate time steps in an implicit method, which is shown to be unconditionally stable.

A related approach was recently developed in [18], where the best rational approximation (BURA) of the implicit time-step operators – including matrix powers – is constructed.

Dealing with the full matrices arising from the discretization of non-local operators is a challenging topic. One can also explore their structure, which makes possible to develop an efficient solution of the corresponding linear problems, see, e.g., [19] and [20].

The aim of the present contribution is to propose an alternative approach, which has the following advances.

•

The analysis for the full discretization error is simple.

•

The spatial convergence order is independent of the power

α

and can ensure a higher-order accuracy.

•

The numerical method is simple and applicable without any extra computations for all exponents

α

^.

To compare which the earlier achievements, note that in the literature, a spatial error of orderh^αappears in the error estimates for any initial conditions, see,e.g., Theorem 3.1 in [15]. We point out that assuming smooth initial condition, this can be changed toh^k, wherekis the polynomial order in the finite element approximation. The smoothness condition will be discussed afterTheorem 2. Also, in the computationally most efficient BURA approach (see [18]), the determination of the coefficients for an arbitrary rational functionr^αneeds efforts.

The main ideas of our works are the following:

•

For the analysis of the spatial discretization error we have used a weak form, which, for the FE discretization is sufficient.

•

The computation algorithm is based on the matrix power–vector product in [21] combined with a conjugated gradient method such that the entire algorithm is composed of sparse matrix–vectors products.

2. Mathematical preliminaries

We investigate the numerical solution of space-fractional diffusion problems. Recall that the (negative) Dirichlet Laplacian operator

−

_∆D

:

L₂(Ω⁾

→

L₂(Ω) is positive and has a compact inverse. The complete orthonormal system of its eigenfunctions and the corresponding eigenvalues will be denoted by{

φ

j

}

j∈Nand{

λ

j

}

j∈N, respectively. Then the fractional Dirichlet Laplacian, which is investigated here, is defined with

Dom (

−

_∆D)^α

= {

u

=

∞

∑

j=₀

u_j

φ

j

:

∞

∑

j=₀

u²_j

· λ

²_j^α

< ∞} ,

⁽

−

_∆D)^αu

=

∞

∑

j=₀

u_j

λ

^α_j

φ

j

.

For more details, we refer to [22].

With this, the space-fractional diffusion equation reads as {

∂

tu(t

,

^x)

= −

(

−

_∆D)^αu(t

,

^{x) x}

∈

_Ω

,

^t

∈

_R⁺

u(0

,

^x)

=

u₀(x) x

∈

_Ω

,

⁽¹⁾

where Ω

⊂

_R^d (d

=

1

,

²

,

3) is a Lipschitz domain. Note that the definition of the Dirichlet Laplacian involves the homogeneous Dirichlet boundary condition.

We use the notion of Sobolev spacesH^k(Ω^{) andH}0^k(Ω) with a non-negative indexk; the corresponding inner product is denoted with

( ·; | · )

kand the notation

∥ · ∥

_kwill be used for the corresponding Sobolev norm. In case ofk

=

0, we usually omit the subscript. To depict clearly the matrix–vector operations in the practical computations, the Euclidian scalar product inR^N will be denoted with (

· , ·

).

In the estimates, the relationr₁≲r₂means thatr₁

≤

cr₂is valid with a positive mesh-independent constantc.

The spatial discretization is performed using a generic finite element spaceV_h^k

⊂

H₀^k(Ω) and the corresponding elliptic projectionP_h

:

H₀^k(Ω⁾

→

V_h^k, which satisfies

(∆DP_h

v | v

h)

=

(∆D

v | v

h)

∀ v

h

∈

V_h^k

.

⁽²⁾

We assume thatV_h^kis chosen so that

∥ v −

P_h

v ∥

_0≲h^k

∥ v ∥

_k

∀ v ∈

H₀^k(Ω⁾

.

⁽³⁾

For the corresponding requirements, we refer to [23, Corollary 1.109]. Using a finite element basis{

ϕ

j

}

j=₁,...,NofV_h^k, one can define the finite element mass matrixM_h

∈

_R^N×N and the stiffness matrixS_h

∈

_R^N×N, the entries of which are given by

[

M_h

]

_i,j

=

(

ϕ

j

| ϕ

i) and

[

S_h

]

_i,j

=

(

ϕ

j

| ϕ

i)₁

=

(

∇ ϕ

j

|∇ ϕ

i)

.

⁽⁴⁾ In the discrete setting, we use the following expansion ofu_h

∈

V_h:

u_h

=

u₁

ϕ

1

+

u₂

ϕ

2

+ · · · +

u_N

ϕ

N (5)

which defines a natural linear bijection

π

h

:

V_h

→

_R^N

,

^with

π

h(u_h)

=

u

=

(u₁

,

^u2

, . . . ,

^uN)^T

.

⁽⁶⁾

For the approximation of the differential operator, we apply the so-called matrix transformation technique. According to this approach, the operation (

−

_∆D)^α

v

is approximated withD^α_h

π

h(

v

h), whereD_h

∈

_R^N×Ndenotes a matrix correspond- ing to

−

_∆DinV_h^k. In case of finite difference discretization, this is straightforward [24], but for finite element methods combined with implicit time-stepping the definition ofD_hneeds a special care.

3. Results

Using the variational principle together with the backward Euler time stepping for

α =

1, the full discretization of(1) can be given as

(uⁿ_h⁺¹

−

uⁿ_h

δ | ϕ

j

)

= −

(

∆^Dun_h⁺¹

| ϕ

j

) j

=

1

,

²

, . . . ,

^N

,

where

δ >

0 is the time step and

uⁿ_h

=

uⁿ₁

ϕ

1

+

uⁿ₂

ϕ

2

+ · · · +

uⁿ_N

ϕ

N

∈

H₀¹(Ω⁾

is the approximation ofu(n

δ, ·

)

:

_Ω

→

_R. According to(4), this can be recasted into the matrix–vector form uⁿ⁺¹

−

uⁿ

δ = −

D_huⁿ⁺¹

,

whereD_h

=

M_h⁻¹S_h. An important observation here, that for using the matrix transformation method, we need to take the power ofD_h instead of the stiffness matrixS_h. With this, the matrix transformation method for a general

α ∈

_R⁺, combined with the time discretization above, can be given as

uⁿ⁺¹

−

uⁿ

δ = −

D^α_huⁿ⁺¹

.

⁽⁷⁾

We perform an error analysis for this approach.

3.1. Spatial discretization

For the error analysis of the spatial discretization, we need the following identities.

Lemma 1. For anya

∈

_R^N and

v

h

∈

V_hwe have

(

π

h⁻¹a

| v

h)

=

(M_ha

, π

h

v

h) (8)

and

(

∇ π

_h^−1a

|∇ v

h)

=

(S_ha

, π

h

v

h)

.

⁽⁹⁾

Proof.

Using the definition of

π

h, and the expansion in(5), we have

(

π

_h^−1a

| v

h)

=

(a₁

ϕ

1

+ · · · +

a_N

ϕ

N

| v

1

ϕ

1

+ · · · + v

N

ϕ

N)

=

M_ha

·

(

v

1

, v

2

, . . . , v

N)^T

=

(M_ha

, π

h

v

h)

,

as stated. The derivation of the second equality can be performed in a similar way. □ For the brevity, we also use the notation

(

−

_∆D,h)^α

= π

_h^−1Dα_h

π

h

,

⁽¹⁰⁾

which gives an approximation of (

−

_∆D)^αonV_h^k.

To obtain a sharp estimate of this term, we use a weak formulation of Balakrishnan’s representation [25] in the case of Hilbert space operators. Henceforth, in the article, we assume that

α ∈

(0

,

^1).

2107

Theorem 1. Let A denote a positive self-adjoint operator on a Hilbert spaceHand

α ∈

(0

,

¹⁾an arbitrary exponent. Then for all u

∈

Dom A and

v ∈

Hwe have the following equality:

(A^αu

| v

)

=

^sin(

πα

⁾

π

∫ ∞ 0

(s^α−1A(sI

+

A)⁻¹u

| v

)

ds

.

⁽¹¹⁾

Henceforth, we analyze the caseA

= −

_∆DandA_h

= −

_∆D,h. The basis of the spatial error estimation is the following statement.

Proposition 1. Using the above notations, we have the following inequality for each eigenfunction

φ

jof

−

_∆D:

⏐

⏐

(((

−

_∆D)^α

−

(

−

_∆D,h)^αP_h)

φ

j

| v

h

)⏐

⏐_≲h^k

λ

^α_j

∥ φ

j

∥

_k

∥ v ∥ .

Proof. Using Balakrishnan’s representation in(11), we have

((A^α

−

A^α_hP_h)

φ

j

| v

h

)

=

^sin(

πα

⁾

π

∫ ∞ 0

s^α−1(

(A(sI

+

A)⁻¹

−

A_h(sI

+

A_h)⁻¹P_h)

φ

j

| v

h

)ds

.

⁽¹²⁾

Using the definition ofP_hin(2), then(9),(8)and(10), we obtain that for an arbitrary

w ∈

H₀^k(Ω) the following identity is valid:

(

^A

w | v

h

) = (

^APh

w | v

h

) = ( ∇

P_h

w |∇ v

h

) = (

^Sh

π

hP_h

w | π

h

v

h

)

= (

^MhD_h

π

hP_h

w | π

h

v

h

) =

(

π

_h^−1Dh

π

hP_h

w | v

h

)

= (

^AhP_h

w | v

h

) .

Using this with

w =

(sI

+

A)⁻¹

φ

j, we can rewrite the scalar product on the right hand side of(12)as follows:

((A(sI

+

A)⁻¹

−

A_h(sI

+

A_h)⁻¹P_h)

φ

j

| v

h

)

=

(

A_h(P_h(sI

+

A)⁻¹

−

(sI

+

A_h)⁻¹P_h)

φ

j

| v

h

)

=

(

(P_h(sI

+

A)⁻¹

−

(sI

+

A_h)⁻¹P_h)

φ

j

|

A_h

v

h

)

.

⁽¹³⁾

Inserting the identity

P_h(sI

+

A)⁻¹

−

(sI

+

A_h)⁻¹P_h

=

(sI

+

A_h)⁻¹

[

(sI

+

A_h)P_h

−

P_h(sI

+

A)

]

(sI

+

A)⁻¹

into the last term of(13), using the self-adjoint property of (sI

+

A_h)⁻¹and(12)again, we obtain that ((A(sI

+

A)⁻¹

−

A_h(sI

+

A_h)⁻¹P_h)

φ

j

| v

h

)

=

(

(sI

+

A_h)⁻¹

[

(sI

+

A_h)P_h

−

P_h(sI

+

A)

]

(sI

+

A)⁻¹

φ

j

|

A_h

v

h

)

=

(

(sI

+

A_h)⁻¹

[

A_hP_h

−

P_hA

]

(sI

+

A)⁻¹

φ

j

|

A_h

v

h

)

=

(

[

A_hP_h

−

P_hA

]

(sI

+

A)⁻¹

φ

j

|

(sI

+

A_h)⁻¹A_h

v

h

)

=

(

[

I

−

P_h

]

A(sI

+

A)⁻¹

φ

j

|

(sI

+

A_h)⁻¹A_h

v

h

)

.

(14)

Applying the approximation property in(3)with(13)and(14)with the Cauchy–Schwarz inequality we finally have that

⏐

⏐

((A(sI

+

A)⁻¹

−

A_h(sI

+

A_h)⁻¹P_h)

φ

j

| v

h

)⏐

⏐

=

⏐

⏐

(

[

I

−

P_h

]

A(sI

+

A)⁻¹

φ

j

|

(sI

+

A_h)⁻¹A_h

v

h

)⏐

⏐

≤ ∥[

I

−

P_h

]

A(sI

+

A)⁻¹

φ

j

∥∥

(sI

+

A_h)⁻¹A_h

v

h

∥ ≤

h^k

∥

A(sI

+

A)⁻¹

φ

j

∥

_k

∥

(sI

+

A_h)⁻¹A_h

v

h

∥

=

h^k

∥ λ

j(sI

+ λ

j)⁻¹

φ

j

∥

_k

∥

(sI

+

A_h)⁻¹A_h

v

h

∥ .

(15)

SinceA_his a positive operator onV_h, we also have

∥

(s

+

A_h)⁻¹A_h

v

h

∥ ≤ ∥ v

h

∥

for everys

≥

0. Inserting this estimate with(15)into the right-hand side of(12), we get

|

(

(A^α

−

A^α_hP_h)

φ

j

| v

h

)

|≤

^sin(

πα

⁾

π

^chk

∫ ∞ 0

s^α−1

λ

j(s

+ λ

j)

∥ φ

j

∥

_k

∥ v

h

∥

ds≲h^k

λ

^αj

∥ φ

j

∥

_k

∥ v

h

∥ ,

as stated in the proposition. □

Using the orthonormal system{

φ

j

}

j∈N, we have the following expansions in(1):

u(0

, ·

)

=

∞

∑

j=₁

u_0j

φ

j

,

^u(t

, ·

)

=

∞

∑

j=₁

u_j

φ

j

,

⁽¹⁶⁾

where{ u_0j}

j∈N

,

{ u_j}

j∈N

⊂

_R, and in the second case, thet-dependence of the coefficientsu_jis not displayed.

Theorem 2. Assume that

∞

∑

j=₀

| λ

^kj^/2+αu_0j

|=

C₀

< ∞

. Then the following estimate holds for a general u

=

u(t

, ·

)in(1):

|

(

(A^α

−

A^α_hP_h)u

| v

h

)

|

_≲C₀h^k

∥ v

h

∥ .

Proof.

The solution of(1)can be given as u(t

, ·

)

=

∞

∑

j=₀

e^−λαj·tu_0,j

φ

j

,

see [26], such that

|

u_j

| < |

u_0j

|

and the assumption imply that

∞

∑

j=0

| λ

^kj^/2+αu_j

| <

^C0 (17)

uniformly for allt

∈

_R⁺.

We also recall a growth condition for theH^k(Ω)-norm of the Dirichlet-Laplacian eigenfunctions:

∥ φ

j

∥

_k≲

λ

_j^k2

.

⁽¹⁸⁾

Note also that by the assumption,u_j

λ

j

→

0. Hence, there is an indexj₀such that for allj

≥

j₀, the following estimate is valid:

(u_j

λ

^αj

)2

≤ |

u_j

λ

^αj

| ≤ λ

_j^k2

|

u_j

λ

^αj

| .

Accordingly, we have

∞

∑

j=₀

(u_j

λ

^α_j)2

< ∞

, and therefore,u

∈

D(A^α) with

A^αu

=

∞

∑

j=₀

u_j

λ

^αj

φ

j

=

∞

∑

j=₀

A^α( u_j

φ

j

)

.

⁽¹⁹⁾

In other words,(19)means for

α =

1 that lim

K→∞A

K

∑

j=₀

u_j

φ

j

=

lim

K→∞

K

∑

j=₀

u_jA

φ

j

=

∞

∑

j=₀

u_jA

φ

j

=

A

∞

∑

j=₀

u_j

φ

j

=

Au

.

In this way, the limit lim

K→∞

K

∑

j=₀

u_j

φ

j

=

u

=

∞

∑

j=₀

u_j

φ

j

is also valid in theH¹(Ω)-norm so that using theH¹(Ω)-orthogonality of the projectionP_hand(16), we finally have P_h

K

∑

j=₀

u_j

φ

j

=

P_hu

=

∞

∑

j=₀

P_h( u_j

φ

j

)

=

lim

K→∞

K

∑

j=₀

P_hu_j

φ

j (20)

Applying the continuous linear operatorA^α_h

:

V_h

→

H₀¹(Ω^{) in}(20), using(19), the estimate inProposition 1, the estimate in(18)and finally the assumption in(17), we have

⏐

⏐

((A^α

−

A^α_hP_h)u

| v

h

)⏐

⏐

=

⏐

⏐

⏐

⏐

⏐

⏐

⎛

⎝(A^α

−

A^α_hP_h)

∞

∑

j=₀

u_j

φ

j

| v

h

⎞

⎠

⏐

⏐

⏐

⏐

⏐

⏐

=

⏐

⏐

⏐

⏐

⏐

⏐

⎛

⎝A^α

∞

∑

j=₀

u_j

φ

j

−

A^α_hP_h

∞

∑

j=₀

u_j

φ

j

| v

h

⎞

⎠

⏐

⏐

⏐

⏐

⏐

⏐

=

⏐

⏐

⏐

⏐

⏐

⏐

⎛

⎝

∞

∑

j=₀

u_jA^α

φ

j

−

∞

∑

j=₀

u_jA^α_hP_h

φ

j

| v

h

⎞

⎠

⏐

⏐

⏐

⏐

⏐

⏐

=

⏐

⏐

⏐

⏐

⏐

⏐

⎛

⎝

∞

∑

j=₀

u_j(A^α

−

A^α_hP_h)

φ

j

| v

h

⎞

⎠

⏐

⏐

⏐

⏐

⏐

⏐

≤

∞

∑

j=₀

⏐

⏐

(u_j(A^α

−

A^α_hP_h)

φ

j

| v

h

)⏐

⏐

≲

∥ v

h

∥

∞

∑

j=₀

⏐

⏐u_j⏐

⏐h^k

λ

^α_j

∥ φ

j

∥

_k≲h^k

∥ v

h

∥

∞

∑

j=₀

⏐

⏐u_j⏐

⏐

λ

^α_j^+k2 ≲h^k

∥ v

h

∥ ,

which proves the statement. □

2109

Remark. If the fractional-order diffusion is observed fromt₀

>

^{0, i.e.u(0}

, ·

)

=

u^∗(t₀

, ·

) for someu^∗, the exponential decay of the coefficientsu_jimplies that∑∞

j=₀(j

+

1)^Ku_0j

< ∞

for any powerK. Therefore, the relation

λ

j

≈

j^d2 immediately gives

∑∞

j=₀

λ

^K_j^u0j

< ∞

, as we have assumed inTheorem 2.

3.2. Time discretization

Now we can prove the convergence of the full discretization. Since for the solutionuof the problem in(1), we have u(t

, ·

)

∈

C^∞(Ω) for every t

∈

(0

,

^T

]

(see Proposition 1 in [9]), we only need the assumption inTheorem 2 for the initial data. Note that regarding the stability, a related result was established in [16] including also second order time discretizations with a possible source term. Here, we also give the convergence rate explicitly.

Theorem 3. Using the assumption inTheorem2for u(0

, ·

)in(1), the full discretization in(7)is convergent of order O(

δ +

h^k), where k is the approximation order of the finite element discretization.

Proof. Rewriting(7)inV_h^kgives the scheme 1

/δ

(

u^j_h⁺¹

−

u^j_h

| v

h

)

+

(

A^α_hu^j_h⁺¹

| v

h

)

=

0

∀ v

h

∈

V_h^k

,

^j

=

0

,

¹

, . . . ,

⁽²¹⁾

whereu^j_his the numerical solution att

=

j

δ

^. □

To analyze this scheme, we use the notationu^j

=

u(j

δ, ·

) withuthe analytic solution of(1)and the elliptic projection P_hintroduced in(2)to obtain the following equality:

1

/δ

(

P_hu^j+1

−

P_hu^j

| v

h

)

+

(

A^α_hP_hu^j+1

| v

h

)

=

1

/δ

(

u^j+1

−

u^j

| v

h

)

+

(

A^αu^j+1

| v

h

)

+

1

/δ

(

(P_hu^j+1

−

u^j+1)

−

(P_hu^j

−

u^j)

| v

h

)

+

(

A^α_hP_hu^j+1

−

A^αu^j+1

| v

h

)

=

(

1

/δ

^(uj+1

−

u^j)

− ∂

tu^j+1

| v

h

)

+

1

/δ

(

(P_hu^j+1

−

u^j+1)

−

(P_hu^j

−

u^j)

| v

h

)

+

(

A^α_hP_hu^j+1

−

A^αu^j+1

| v

h

)

:=

(z^j

| v

h)

,

(22)

wherez^jcan be recognized as a consistency error, which we estimate termwise.

Obviously, using the general mean value theorem, and the smoothness of the analytic solution, we obtain

∥

1

/δ

^(uj+1

−

u^j)

− ∂

tu^j+1

∥ =







 1

/δ

∫ _δ

0

∂

tu(j

δ +

s

, ·

)

− ∂

tu((j

+

1)

δ, ·

) ds









≤

1

/δ

∫ _δ

0

∥ ∂

tu(j

δ +

s

, ·

)

− ∂

tu((j

+

1)

δ, ·

)

∥

ds

≤

1

/δ

∫ _δ

0

δ

^max

ξ∈[₀,δ]

∥ ∂

ttu(j

δ + ξ, ·

)

∥

ds≲

δ.

Similarly, the general mean value theorem and the approximation property in(3)imply

∥

1

/δ

^((Phu^j+1

−

u^j+1)

−

(P_hu^j

−

u^j))

∥ ≤

ch^k

∥ ∂

tu

∥

_L∞([_jk,^(j+_1)k];_H^k₍Ω⁾⁾

.

Combining these estimates withTheorem 2, we get

|

( z^j

| v

h

)

|≤

(O(

δ

⁾

+

O(h^k+1))

∥ v

h

∥ .

⁽²³⁾

Using the notationy^j_h

=

P_hu^j

−

u^j_hin the difference of(21)and(22), we obtain that 1

/δ

(

y^j_h⁺¹

−

y^j_h

| v

h

)

+

(

A^α_hy^j_h⁺¹

| v

h

)

=

( z^j

| v

h

)

∀ v

h

∈

V_h

.

Taking here

v

h

=

y^j_h⁺¹and rearranging the equality, we have

∥

y^j_h⁺¹

∥

²

+ δ

(

A^α_hy^j_h⁺¹

|

y^j_h⁺¹)

=

(

y^j_h

+ δ

^zj

|

y^j_h⁺¹)

,

and therefore, using the positivity ofA^α_h and(23), we get

∥

y^j_h⁺¹

∥

²

≤ ∥

y^j_h

∥∥

y^j_h⁺¹

∥ + δ

( z^j

|

y^j_h⁺¹)

≤ ∥

y^j_h

∥∥

y^j_h⁺¹

∥ + δ ·

O(

δ +

h^k)

∥

y^j_h⁺¹

∥ ,

which results in the following estimate:

∥

y^j_h⁺¹

∥ ≤ ∥

y^j_h

∥ + δ

^(O(

δ +

h^k))

.

A consecutive application of this inequality gives max

0≤_j≤_M

∥

y^j_h

∥ ≤ ∥

y⁰_h

∥ +

T

·

O(

δ +

h^k)

,

which, together with(3)can be used to get the final error estimator max

0≤_j≤_M

∥

u^j

−

u^j_h

∥ ≤

max

0≤_j≤_M

(

∥

u^j

−

P_hu^j

∥ + ∥

P_hu^j

−

u^j_h

∥

)

=

max

0≤_j≤_M

(

∥

u^j

−

P_hu^j

∥ + ∥

y^j_h

∥

)

=

O(

δ +

h^k)

,

as stated in the theorem.

4. Numerical method

For the numerical solution of problem(1), we need to use an efficient method avoiding the direct computation of the

α

-th power of the matrices. Such an approach was first proposed in [21], which can be applied immediately only for explicit time stepping. Here we propose an algorithm to approximate the implicit time steppingu^j+1

=

(I

+ δ

^Aα_h^)−1uj, i.e.uⁿ

=

(I

+ δ

^Aα_h^)−nu0without computing matrix powers or solving linear systems with a dense matrix. We also assume here thatA_his symmetric, which will be satisfied for the finite element space in the numerical experiments.

In a general situation, it is also satisfied if the finite element basis functions are translations of each other. In practice, we can use this if the domain is approximated using a square grid.

In a general situation, for arbitrary finite elements, anL²-orthogonal basis impliesM_h

=

I_h such that A_h

=

I_h⁻¹S_h becomes symmetric. This, however, leads to a full matrixA_h, which will slow down our algorithm.

For this, we will combine the algorithm in [21] with a conjugate gradient method.

4.1. The algorithm

The proposed method consists of the following steps.

(i) Following the method in [21], we compute the k₁ smallest and k₂ largest eigenvalues. In an increasing order with_k

¯ =

k₁

+

k₂these and the corresponding eigenvectors are denoted by

λ

1

, λ

2

, . . . , λ

k^¯ ofA_h andx₁

,

^x2

, . . . ,

^x¯k, respectively.

(ii) LetX

= [

x₁

,

^x2

, . . . ,

^x_k^¯

] ∈

_R^nׯkdenote the matrix composed of these eigenvectors, andQ_k¯

=

XX^T the orthogonal projection matrix to the subspace span

{

x₁

,

^x2

, . . . ,

^xk^¯

}

. In this case,I

−

Q_k¯is the orthogonal projection matrix to the complementary subspace.

(iii) The problem is divided then into two parts:

(I

+ δ

^Aαh)⁻ⁿu^j

=

(I

+ δ

^Aαh)⁻ⁿQ_k¯u^j

+

(I

+ δ

^Aαh)⁻ⁿ(I

−

Q¯_k)u^j (24) (a) We can directly compute the first part, since we already know the corresponding eigenvalues:

(I

+ δ

^Aαh)⁻ⁿQ¯_ku^j

=

XΛ^nXTuj

,

whereΛdenotes the diagonal matrix consisting of the elements ₁₊¹_δλα

1

,

₁+¹δλ^α2

, . . . ,

₁+¹δλ^α_k¯.

(b) To approximate (I

+ δ

^Aα_h^)−n(I

−

Q_k¯)u^j, we apply a conjugate gradient methodCG, see,e.g.[], such that CGw

≈

(I

+

A^α_h)⁻¹w

.

In the steps of the conjugate gradient algorithm, we use an approximation of the matrix–vector products (I

+ δ

^Aα_h^)w without computingA^α_h. This is performed using the Taylor series approach

(2A_h

σ

^(A)

)α

w

=

∞

∑

n=₀

(

α

n

) (2A_h

σ

^(A)

−

I

)n

w

≈

K

∑

n=₀

(

α

n

) (2A_h

σ

^(A)

−

I

)n

w

:=

( 2

σ

^(A)

)α

T(A_h

, α

^)w

,

⁽²⁵⁾ where

σ

(A) denotes the spectral radius ofA.

Remark. The conjugate gradient algorithm is suitable here as we have a symmetric positive definite matrixI

+

A^α_h. The stopping criterion (or tolerance) for this procedure is given by discussing the numerical experiments.

An important parameter in the approximation in(25)is the parameterK. An estimation of this is given in Section2.3 in [21], which motivated our choice in Section4.3.

The operations of the conjugate gradient method are invariant to the subspace ran(I

−

Q¯_k), thus the Taylor-series method will converge quickly in every time step. Also, in the Taylor series approach, we use only sparse matrix–vector products such that beyond the eigenvalue approximation, the entire algorithm involves only these very quick operations.

4.2. Error analysis of the algorithm

Recall that inTheorem 3, we estimated the difference between the analytic solution and the numerical solution based on the implicit Euler time stepping. In the practice, however, according to (iii)(b) in the above algorithm, we do not apply directly the implicit time steps. In concrete terms, we computeCGⁿw⁰instead of (I

+ δ

^Aαh)⁻ⁿw⁰, wherew⁰

=

(I

−

Q_k¯)u⁰.

2111