DeepAI

# Galerkin Type Methods for Semilinear Time-Fractional Diffusion Problems

We derive optimal L^2-error estimates for semilinear time-fractional subdiffusion problems involving Caputo derivatives in time of order α∈ (0,1), for cases with smooth and nonsmooth initial data. A general framework is introduced allowing a unified error analysis of Galerkin type space approximation methods. The analysis is based on a semigroup type approach and exploits the properties of the inverse of the associated elliptic operator. Completely discrete schemes are analyzed in the same framework using a backward Euler convolution quadrature method in time. Numerical examples including conforming, nonconforming and mixed finite element (FE) methods are presented to illustrate the theoretical results.

12/04/2020

### Numerical solution of a time-fractional nonlinear Rayleigh-Stokes problem

We study a semilinear fractional-in-time Rayleigh-Stokes problem for a g...
07/31/2020

### Numerical Analysis of Backward Subdiffusion Problems

The aim of this paper is to develop and analyze numerical schemes for ap...
09/09/2021

### Fractional Crank-Nicolson-Galerkin finite element methods for nonlinear time fractional parabolic problems with time delay

A linearized numerical scheme is proposed to solve the nonlinear time fr...
09/28/2019

### Complete monotonicity-preserving numerical methods for time fractional ODEs

The time fractional ODEs are equivalent to the convolutional Volterra in...
10/06/2020

### Error Analysis of Finite Element Approximations of Diffusion Coefficient Identification for Elliptic and Parabolic Problems

In this work, we present a novel error analysis for recovering a spatial...
06/15/2019

### Incomplete Iterative Solution of the Subdiffusion Problem

In this work, we develop an efficient incomplete iterative scheme for th...
01/17/2023

### Error estimates for completely discrete FEM in energy-type and weaker norms

The paper presents error estimates within a unified abstract framework f...

## 1 Introduction

The purpose of this paper is to discuss some aspects of the numerical solution of the semilinear time-fractional initial boundary value problem

 C∂αtu+Lu=f(x,t,u) in Ω×(0,T0],u(x,0)=u0(x) in Ω, (1.1)

subject to a homogeneous Dirichlet boundary condition, where () is a bounded convex polyhedral domain with a boundary and is a fixed time. Here is a given initial data and is a smooth function of its arguments satisfying

 supx∈Ω,t∈(0,T0)(|∂tf(x,t,u)|+|∂uf(x,t,u)|)≤L∀u∈R. (1.2)

The operator is defined by , where is a symmetric and uniformly positive definite in matrix, and is nonnegative. The coefficients and are assumed to be sufficiently smooth on . The operator is the Caputo fractional derivative in time of order defined by

 C∂αtφ(t)=1Γ(1−α)∫t0(t−s)−α∂sφ(s)ds,0<α<1, (1.3)

where and denotes the usual Gamma function. As , converges to , and thus, problem (1.1) reduces to the standard semilinear parabolic problem [34].

Let denote the inner product in with induced norm . Since is convex, the solution of the elliptic problem in , with on and , belongs to . With , recall that the operator is selfadjoint, positive definite and has a compact inverse. Let

denotes the eigenvalues and eigenfunctions of

with an orthonormal basis in . By spectral method, the fractional powers of are defined by

 Lνv=∞∑j=1λνj(v,φj)φj,ν>0,

with domains . Note that

form a Hilbert scale of interpolation spaces and

with continuous and compact embeddings for .

The regularity of the solution in (1.1) plays a key role in our error analysis. For initial data , , problem (1.1) has a unique solution satisfying [1, Theorem 3.1]:

 u∈Cαν([0,T0];L2(Ω))∩C([0,T0];D(Lν))∩C((0,T0];D(L)), (1.4)
 C∂αtu∈C((0,T0];L2(Ω)), (1.5)
 ∂tu(t)∈L2(Ω)and∥∂tu(t)∥≤ctαν−1,t∈(0,T0]. (1.6)

The results show that the solution of the semilinear problem (1.1) enjoys (to some extent) smoothing properties analogous to those of the homogeneous linear problem. For , it is shown that ([1, Theorem 3.2])

 u∈C([0,T0];L2(Ω))∩Lγ(0,T0;D(L)),γ<1/α. (1.7)

Note that the first time derivative of is not smooth enough in space even in the case of a smooth initial data. This actually causes a major difficulty in deriving optimal error estimates based on standard techniques, such as the energy method.

The numerical approximation of fractional differential equations has received considerable attention over the last two decades. For linear time-fractional equations, a vast literature is now available. See the short list [26, 27, 11, 10, 16, 18] on problems with nonsmooth data and [12] for a concise overview and recent developments. In contrast, numerical studies on nonlinear time-fractional evolution problems are rather limited. In [21], a linearized -Galerkin FEM was proposed to solve a nonlinear time-fractional Schrödinger equation. In [20], -type schemes have been analyzed for approximating the solution of (1.1). The error estimates in [21] and [20] are derived under high regularity assumptions on the exact solution, so the limited smoothing property of the model (1.1) was not taken into consideration. In [13], the numerical solution of (1.1) was investigated assuming that the nonlinearity is uniformly Lipschitz in and the initial data . Error estimates are established for linearized time-stepping schemes based on the -method and a convolution quadrature generated by the backward Euler difference formula. In the recent paper [1], we derived error estimates for the same problem with initial data , . The new estimates extend known results obtained for the standard semilinear parabolic problem [14]. For other types of time-fractional problems, one may refer to [5] for fractional diffusion-wave equations and to [28] for an integro-differential equation.

In this paper, we approximate the solution of the semilinear problem (1.1) by general Galerkin type approximation methods in space and a convolution quadrature in time. Our aim is to develop a unified error analysis with optimal error estimates with respect to the data regularity. We shall follow a semigroup type approach and make use of the inverse of the associated elliptic operator [3]. The current study extends the recent work [16] dealing with the homogeneous linear problem, which relied on the energy technique. Our analysis includes conforming, nonconforming and mixed FEMs, and the results are applicable to nonlinear multi-term diffusion problems. It is worth noting that most of our results hold in the limiting case , i.e., our study also generalizes the work [3]. Particularly interesting are the estimates derived for the mixed FEM, which are new and have not been established earlier.

The paper is organized as follows. In section 2, a general setting of the problem is introduced and preliminary error estimates are derived, which require regularity properties analogous to those of the homogeneous linear problem. In section 3, an alternative error estimation is proposed without a priori regularity assumptions on the exact solution. Time-stepping schemes based on a backward Euler convolution quadrature method are analyzed in section 4. Applications are presented in section 5. The mixed form of problem (1.1) is considered in section 6 and related convergence rates are obtained. Finally, numerical results are provided to validate the theoretical findings.

Throughout the paper, we denote by a constant which may vary at different occurrences, but is always independent of the mesh size and the time step size . We shall also use the abbreviation and for and , respectively.

## 2 General setting and preliminary estimates

Set . Then, is compact, selfadjoint and positive definite. In terms of , we may write (1.1) as

 TC∂αtu+u=Tf(u),t>0,u(0)=u0. (2.1)

For the purpose of approximating the solution of this problem, let be a family of finite-dimensional spaces that depends on , . We assume that we are given a corresponding family of linear operators which approximate . Then consider the semidiscrete problem: find for such that

 TCh∂αtuh+uh=Thf(uh),t>0,uh(0)=u0h∈Vh, (2.2)

where is a suitably chosen approximation of . In our analysis, we shall make the following assumptions:
(i) is selfadjoint, positive semidefinite on and positive definite on .
(ii) , where is the orthogonal -projection onto .
(iii) For some constants and , there holds

 ∥Thf−Tf∥≤chγ∥f∥∀f∈L2(Ω). (2.3)

Since exists on , (2.2) may be solved uniquely for . The following diagram displays the different links between the operators under consideration:

In the diagram, the operator is defined by . It is the analogue of the Ritz elliptic projection in the context of Galerkin FE methods. Note that , and in view of (2.3), satisfies

 ∥Rhv−v∥=∥ThLv−TLv∥≤chγ∥Lv∥∀v∈D(L). (2.4)

Further, by the definition of , we see that .

Examples of family with the above properties are exhibited by the standard Galerkin FE and spectral methods in the case , and by other nonconforming Galerkin methods in the case . The mixed FE method applied to (1.1) is a typical example which has the above properties and will be considered in this study.

By our assumptions on , the operator satisfies

 ∥(z−αI+Th)−1∥≤M|z|α∀z∈Σθ, (2.5)

where is the sector with being fixed and depends on . In (2.5), and in the sequel, we keep the same notation to denote the operator norm from . Using that

 (z−αI+Th)−1Th=I−z−α(z−αI+Th)−1, (2.6)

we obtain

 ∥(z−αI+Th)−1Th∥≤1+M∀z∈Σθ. (2.7)

Note that (2.5) and (2.7) hold for . By means of the Laplace transform, the solution of problem (2.2) is represented by

 uh(t)=Eh(t)u0h+∫t0¯Eh(t−s)f(uh(s))ds,t>0. (2.8)

The operators and are defined by

 Eh(t)=12πi∫Γθ,δeztz−1Kh(z)dz and ¯Eh(t)=12πi∫Γθ,δeztz−αKh(z)dz, (2.9)

respectively, where . The contour with and , is oriented with an increasing imaginary part. Similarly, the solution of problem (2.1) is given by

 u(t)=E(t)u0+∫t0¯E(t−s)f(u(s))ds,t>0, (2.10)

where the operators and are respectively defined in terms of as in (2.9). Standard arguments show that (see for instance [25])

 ∥E(t)v∥+∥Eh(t)v∥+t1−α(∥¯E(t)v∥+∥¯Eh(t)v∥)≤c∥v∥∀v∈L2(Ω). (2.11)

Now let denote the error at time . Define the intermediate solution , , by

 TCh∂αtvh+vh=Thf(u),t>0,vh(0)=u0h. (2.12)

Then, by splitting the error , and subtracting (2.12) from (2.1), we find that satisfies

 TCh∂αtξ(t)+ξ(t)=(Th−T)(f(u)−C∂αtu)(t),t>0. (2.13)

With

 ρ(t):=(Th−T)(f(u)−C∂αtu)(t), (2.14)

we thus obtain

 TCh∂αtξ(t)+ξ(t)=ρ(t),t>0,ξ(0)∈L2(Ω). (2.15)

Before proving the main result of this section, we recall the following lemma which generalizes the classical Gronwall’s inequality, see [6].

###### Lemma 2.1.

Assume that is a nonnegative function in which satisfies

 y(t)≤g(t)+β∫t0(t−s)−αy(s)ds for t∈(0,T0],

where , , and . Then there exists a constant such that

 y(t)≤g(t)+CT0∫t0(t−s)−αg(s)ds for t∈(0,T0].

Note that, by using (2.10), (2.11) and the inequality

 ∥f(u(t))∥≤∥f(u(t))−f(0)∥+∥f(0)∥≤L∥u∥+∥f(0)∥,t≥0, (2.16)

Lemma 2.1 implies that for with . Now we are ready to prove an error estimate for problem (2.1). Here is given by (2.14) and .

###### Lemma 2.2.

Let and be the solutions of (2.1) and (2.2), respectively. Assume that . Then, for ,

 ∥e(t)∥≤c(G(t)+∫t0(t−s)α−1G(s)ds), (2.17)

where

 G(t)=t−1sups≤t(∥~ρ(s)∥+s∥ρ(s)∥+s2∥ρt(s)∥), (2.18)

and is independent of .

###### Proof.

First we derive a bound for the difference between and the intermediate solution . Since , an application of Lemma 3.5 in [16] to (2.1) and (2.12) yields , where is given by (2.18). In view of the splitting , it suffices to estimate . Note that satisfies

 TCh∂αtη(t)+η(t)=Th(f(uh)−f(u)),t>0,η(0)=0.

Hence, by Duhamel’s principle,

 η(t)=∫t0¯Eh(t−s)(f(uh(s))−f(u(s)))ds,t>0.

Using the property of in (2.11) and condition (1.2), we see that

 ∥η(t)∥≤cL∫t0(t−s)α−1∥uh(s)−u(s)∥ds,t>0,

and thus

 ∥e(t)∥≤∥ξ(t)∥+c∫t0(t−s)α−1∥e(s)∥ds,t>0.

An application of Lemma 2.1 yields (2.17), which completes the proof. ∎

Clearly the error estimate in Lemma 2.2 is meaningful provided that . Recalling that , we have by (2.3), . Hence, to achieve a order of convergence, we need to assume that for . It turns out that, without additional conditions on initial data and nonlinearity, this property, which holds in the linear case, does not generalize to the semilinear problem. This remark equally applies to the semilinear parabalic problem, see the discussion in [31, pp. 228].

## 3 Error estimates without regularity assumptions

We shall present below an alternative derivation of the error bound without a priori regularity assumptions on the exact solution. To do so, we first introduce the operator

 Sh(z)=(z−αI+Th)−1Th−(z−αI+T)−1T.

Then satisfies the following property.

###### Lemma 3.1.

There holds

 ∥Sh(z)v∥≤chγ|z|α(1−ν)∥Lνv∥∀z∈Σθ,ν∈[0,1]. (3.1)
###### Proof.

Using the identity (2.6), we verify that

 Sh(z) = z−α(z−αI+T)−1[(z−αI+Th)−(z−αI+T)](z−αI+Th)−1 = z−α(z−αI+T)−1(Th−T)(z−αI+Th)−1.

Then, by (2.5) and (2.3),

 ∥Sh(z)v∥≤c∥(Th−T)(z−αI+Th)−1v∥≤chγ|z|α∥v∥.

This shows (3.1) for . For , i.e., , we have . Then, by (2.3) and (2.7), we get

 ∥Sh(z)v∥≤c∥(Th−T)(z−αI+T)−1TLv∥≤chγ∥Lv∥.

The desired estimate (2.3) follows now by interpolation. ∎

We further introduce the following operators: and . Then, by Lemma 3.1,

 ∥¯Fh(t)v∥≤chγ∫Γθ,1/teRe(z)t|dz|∥v∥≤ct−1hγ∥v∥. (3.2)

Similarly, based on (3.1), the following estimate

 ∥Fh(t)v∥≤ct−α(1−ν)hγ∥Lνv∥,ν=0,1, (3.3)

holds for . Now we are ready to prove a nonsmooth data error estimate. Here and the throughout the paper, if and otherwise.

###### Theorem 3.1.

Let , . Let and be the solutions defined by (2.10) and (2.8), respectively, with . Then there is a constant , where , such that

 ∥uh(t)−u(t)∥≤chγℓh(ν)t−α(1−ν),t∈(0,T0]. (3.4)
###### Proof.

Recall that . From (2.8) and (2.10), we get after rearrangements

 e(t)=Fh(t)u0+∫t0¯Eh(t−s)[f(uh(s))−f(u(s))]ds+∫t0¯Fh(t−s)f(u(s))ds. (3.5)

The last term in (3.5) can be written as where

 I=∫t0¯Fh(t−s)(f(u(s))−f(u(t)))dsandII=(∫t0¯Fh(t−s)ds)f(u(t)).

For , we use (3.2), (1.2) and the property to get

 ∥I∥≤chγ∫t0(t−s)−1(t−s)ανds≤chγtαν.

To estimate , we introduce the operator Then and for all since . Hence, as is bounded, see (2.16). Now, using the properties of and in (2.11) and (3.3), respectively, we obtain

 ∥e(t)∥≤ct−α(1−ν)hγ∥Lνu0∥+c∫t0(t−s)α−1∥e(s)∥ds+chγ+chγtαν. (3.6)

The desired estimate follows now by applying Lemma 2.1. To establish the estimate for , we follow the same arguments presented in the proof of Theorem 4.4 in [1]. ∎

As an immediate application of Theorem 3.1, consider the standard conforming Galerkin FEM with consists of piecewise linear functions on a shape-regular triangulation with a mesh parameter . Let be the solution operator of the discrete problem:

 a(Thf,v)=(f,v)∀v∈Vh, (3.7)

where is the bilinear form associated with the elliptic operator . Then, is selfadjoint, positive semidefinite on and positive definite on , see [3], and satisfies (2.3) with . Thus, by Theorem 3.1,

 ∥uh(t)−u(t)∥≤ch2ℓh(ν)t−α(1−ν),t>0, (3.8)

for . This improves the following estimate

 ∥uh(t)−u(t)∥≤ch2(t−α(1−ν)+max(0,ln(tα(1−ν)/h2))),t>0, (3.9)

established in [1]. We notice that the logarithmic factor is also present in the parabolic case when , see [14, Theorem 1.1].

As a second example, we show that the present semidiscrete error analysis extends to the following multi-term time-fractional diffusion problem:

 P(∂t)u+Lu=f(u) in Ω×(0,T0],u(0)=u0 in Ω,u=0 on ∂Ω×(0,T0], (3.10)

where the multi-term differential operator is defined by with being the orders of the fractional Caputo derivatives, and , . This model was derived to improve the modeling accuracy of the single-term model (1.1) for describing anomalous diffusion. An inspection of the proof of the Theorem 3.1 reveals that its main arguments are based on the bounds derived for the operators , and . Following [11], one can verify that these operators satisfy the same bounds as in the single-term case. This readily implies that the estimate (3.4) remains valid for the multi-term diffusion problem (3.10).

###### Remark 3.1.

In the parabolic case, a singularity in time appears which has the same form as in (3.2). Hence, the estimate (3.4) remains valid when .

###### Remark 3.2.

For smooth initial data , the estimate (3.4) still holds for the choice . Indeed, we have

 Eh(t)Rhu0−E(t)u0=Eh(t)(Rhu0−u0)+(Eh(t)u0−E(t)u0).

By the stability of the operator ,

 ∥Eh(t)(Rhu0−u0)∥≤c∥Rhu0−u0∥≤chγ∥Lu0∥.

Then we reach our conclusion by following the arguments in the proof of Theorem 3.1.

## 4 Fully discrete schemes

This section is devoted to the analysis of a fully discrete scheme for problem (2.2) based on a convolution quadrature (CQ) generated by the backward Euler method, using the framework developed in [25, 5]. Divide the time interval into equal subintervals with a time step size , and let . The convolution quadrature [23] refers to an approximation of any function of the form as

 (k∗φ)(tn):=∫tn0k(tn−s)φ(s)ds≈n∑j=0βn−j(τ)φ(tj),

where the weights are computed from the Laplace transform of rather than the kernel . With being time differentiation, define as the operator of (distributional) convolution with the kernel : for a function with suitable smoothness. Then a convolution quadrature will approximate by a discrete convolution at as where the quadrature weights are determined by the generating power series with being a rational function, chosen as the quotient of the generating polynomials of a stable and consistent linear multistep method. For the backward Euler method, .

An important property of the convolution quadrature is that it maintains some relations of the continuous convolution. For instance, the associativity of convolution is valid for the convolution quadrature [5] such as

 K2(∂τ)K1(∂τ)=K2K1(∂τ) and K2(∂τ)(k1∗φ)=(K2(∂τ)k1)∗φ. (4.1)

In the following lemma, we state an interesting result on the error of the convolution quadrature [24, Theorem 5.2].

###### Lemma 4.1.

Let be analytic in the sector and such that

 ∥G(z)∥≤c|z|−μ∀z∈Σθ,

for some real and . Then, for , the convolution quadrature based on the backward Euler method satisfies

 (4.2)

Upon using the relation between the Riemann-Liouville derivative denoted by and the Caputo derivative , the semidiscrete scheme (2.2) can be rewritten as

 Th∂αt(uh−u0h)+uh=Thf(uh),t>0,uh(0)=u0h. (4.3)

Thus, the proposed backward Euler CQ scheme is to seek , , such that

 Th∂ατ(Unh−U0h)+Unh=Thf(Unh),n≥1,U0h=u0h. (4.4)

### 4.1 The linear case

We begin by investigating the time discretization of the inhomogeneous linear problem

 TC∂αtu+u=Tf(t),t>0,u(0)=u0, (4.5)

with a semidiscrete solution satisfying

 TCh∂αtuh+uh=Thf(t),t>0,uh(0)=Phu0. (4.6)

The fully discrete solution is defined by

 Th∂ατ(Unh−U0h)+Unh=Thf(tn),n≥1,U0h=Phu0. (4.7)

Then we establish the following result.

###### Theorem 4.1.

Let , . Let and be the solutions of and , respectively, with . Then, for ,

 ∥Unh−u(tn)∥≤c(τtαν−1n+hγt−α(1−ν)n)∥Lνu0∥. (4.8)
###### Proof.

We first notice that , which follows from Theorem 3.1 in the case . In order to estimate , apply and to (4.6) and (4.7), respectively, use the associativity of convolution and the property , to deduce that

 Unh−uh(tn)=−(G(∂τ)−G(∂t))u0, (4.9)

where . We recall that, by (2.5), Then, Lemma 4.1 (with , ) and the -stability of yield

 ∥Unh−uh(tn)∥≤cτt−1n∥u0∥. (4.10)

For , consider first the choice . Recalling that , we use the identity