Galerkin Type Methods for Semilinear Time-Fractional Diffusion Problems

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} \partial_{t}^{α} u + L u = f (x, t, u) in Ω \times (0, T_{0}], u (x, 0) = u_{0} (x) in Ω,

(1.1)

subject to a homogeneous Dirichlet boundary condition, where $Ω \subset R^{d}$ ( $d \geq 2$ ) is a bounded convex polyhedral domain with a boundary $\partial Ω$ and $T_{0} > 0$ is a fixed time. Here $u_{0}$ is a given initial data and $f$ is a smooth function of its arguments satisfying

sup x \in Ω, t \in (0, T_{0}) (| \partial_{t} f (x, t, u) | + | \partial_{u} f (x, t, u) |) \leq L \forall u \in R .

(1.2)

The operator $L$ is defined by $L u = - div [A (x) \nabla u] + κ (x) u$ , where $A (x) = [a_{i j} (x)]$ is a $d \times d$ symmetric and uniformly positive definite in $¯ Ω$ matrix, and $κ \in L^{\infty} (Ω)$ is nonnegative. The coefficients $a_{i j}$ and $κ$ are assumed to be sufficiently smooth on $¯ Ω$ . The operator $^{C} \partial_{t}^{α}$ is the Caputo fractional derivative in time of order $α \in (0, 1)$ defined by

^{C} \partial_{t}^{α} φ (t) = \frac{1}{Γ (1 - α)} \int_{0}^{t} (t - s)^{- α} \partial_{s} φ (s) d s, 0 < α < 1,

(1.3)

where $\partial_{s} φ = \partial φ / \partial s$ and $Γ (\cdot)$ denotes the usual Gamma function. As $α \to 1^{-}$ , $^{C} \partial_{t}^{α}$ converges to $\partial_{t}$ , and thus, problem (1.1) reduces to the standard semilinear parabolic problem [34].

Let $(\cdot, \cdot)$ denote the inner product in $L^{2} (Ω)$ with induced norm $∥ \cdot ∥$ . Since $Ω$ is convex, the solution of the elliptic problem $L u = f$ in $Ω$ , with $u = 0$ on $\partial Ω$ and $f \in L^{2} (Ω)$ , belongs to $H^{2} (Ω)$ . With $D (L) = H^{2} (Ω) \cap H_{0}^{1} (Ω)$ , recall that the operator $L : D (L) \to L^{2} (Ω)$ is selfadjoint, positive definite and has a compact inverse. Let ${λ_{j}, φ_{j}}_{j = 1}^{\infty}$

denotes the eigenvalues and eigenfunctions of

L

with

{φ_{j}}_{j = 1}^{\infty}

an orthonormal basis in

L^{2} (Ω)

. By spectral method, the fractional powers of

L

are defined by

L^{ν} v = \infty \sum j = 1 λ_{j}^{ν} (v, φ_{j}) φ_{j}, ν > 0,

with domains $D (L^{ν}) = {v \in L^{2} (Ω) : ∥ L^{ν} v ∥ < \infty}$ . Note that ${D (L^{ν})}$

form a Hilbert scale of interpolation spaces and

D (L) \subset D (L^{ν}) \subset D (L^{β}) \subset D (L^{0}) = L^{2} (Ω)

with continuous and compact embeddings for

0 < β < ν < 1

The regularity of the solution in (1.1) plays a key role in our error analysis. For initial data $u_{0} \in D (L^{ν})$ , $ν \in (0, 1]$ , problem (1.1) has a unique solution $u$ satisfying [1, Theorem 3.1]:

u \in C^{α ν} ([0, T_{0}]; L^{2} (Ω)) \cap C ([0, T_{0}]; D (L^{ν})) \cap C ((0, T_{0}]; D (L)),

(1.4)

^{C} \partial_{t}^{α} u \in C ((0, T_{0}]; L^{2} (Ω)),

(1.5)

\partial_{t} u (t) \in L^{2} (Ω) and ∥ \partial_{t} u (t) ∥ \leq c t^{α ν - 1}, t \in (0, T_{0}] .

(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 $u_{0} \in L^{2} (Ω)$ , it is shown that ([1, Theorem 3.2])

u \in C ([0, T_{0}]; L^{2} (Ω)) \cap L^{γ} (0, T_{0}; D (L)), γ < 1 / α .

(1.7)

Note that the first time derivative of $u$ 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 $L^{1}$ -Galerkin FEM was proposed to solve a nonlinear time-fractional Schrödinger equation. In [20], $L^{1}$ -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 $f$ is uniformly Lipschitz in $u$ and the initial data $u_{0} \in D (L)$ . Error estimates are established for linearized time-stepping schemes based on the $L^{1}$ -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 $u_{0} \in D (L^{ν})$ , $ν \in (0, 1]$ . 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 $α = 1$ , 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 $c$ a constant which may vary at different occurrences, but is always independent of the mesh size $h$ and the time step size $τ$ . We shall also use the abbreviation $f (u)$ and $f (t)$ for $f (x, t, u)$ and $f (x, t)$ , respectively.

2 General setting and preliminary estimates

Set $T = L^{- 1}$ . Then, $T : L^{2} (Ω) \to D (L)$ is compact, selfadjoint and positive definite. In terms of $T$ , we may write (1.1) as

T^{C} \partial_{t}^{α} u + u = T f (u), t > 0, u (0) = u_{0} .

(2.1)

For the purpose of approximating the solution of this problem, let $V_{h} \subset L^{2} (Ω)$ be a family of finite-dimensional spaces that depends on $h$ , $0 < h < 1$ . We assume that we are given a corresponding family of linear operators $T_{h} : L^{2} (Ω) \to V_{h}$ which approximate $T$ . Then consider the semidiscrete problem: find $u_{h} (t) \in V_{h}$ for $t \geq 0$ such that

T_{h}^{C} \partial_{t}^{α} u_{h} + u_{h} = T_{h} f (u_{h}), t > 0, u_{h} (0) = u_{0 h} \in V_{h},

(2.2)

where $u_{0 h}$ is a suitably chosen approximation of $u_{0}$ . In our analysis, we shall make the following assumptions:
(i) $T_{h}$ is selfadjoint, positive semidefinite on $L^{2} (Ω)$ and positive definite on $V_{h}$ .
(ii) $T_{h} P_{h} = T_{h}$ , where $P_{h} : L^{2} (Ω) \to V_{h}$ is the orthogonal $L^{2}$ -projection onto $V_{h}$ .
(iii) For some constants $γ > 0$ and $c > 0$ , there holds

∥ T_{h} f - T f ∥ \leq c h^{γ} ∥ f ∥ \forall f \in L^{2} (Ω) .

(2.3)

Since $T_{h}^{- 1}$ exists on $V_{h}$ , (2.2) may be solved uniquely for $t > 0$ . The following diagram displays the different links between the operators under consideration:

In the diagram, the operator $R_{h} : D (L) \to V_{h}$ is defined by $R_{h} = T_{h} L$ . It is the analogue of the Ritz elliptic projection in the context of Galerkin FE methods. Note that $R_{h} T = T_{h}$ , and in view of (2.3), $R_{h}$ satisfies

∥ R_{h} v - v ∥ = ∥ T_{h} L v - T L v ∥ \leq c h^{γ} ∥ L v ∥ \forall v \in D (L) .

(2.4)

Further, by the definition of $P_{h}$ , we see that $∥ P_{h} v - v ∥ \leq ∥ R_{h} v - v ∥ \forall v \in D (L)$ .

Examples of family ${T_{h}}$ with the above properties are exhibited by the standard Galerkin FE and spectral methods in the case $V_{h} \subset H_{0}^{1} (Ω)$ , and by other nonconforming Galerkin methods in the case $V_{h} \subset̸ H_{0}^{1} (Ω)$ . 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 $T_{h}$ , the operator $(z^{- α} I + T_{h})^{- 1} : L^{2} (Ω) \to L^{2} (Ω)$ satisfies

∥ (z^{- α} I + T_{h})^{- 1} ∥ \leq M | z |^{α} \forall z \in Σ_{θ},

(2.5)

where $Σ_{θ}$ is the sector $Σ_{θ} = {z \in C, z \neq 0, | arg z | < θ}$ with $θ \in (π / 2, π)$ being fixed and $M$ depends on $θ$ . In (2.5), and in the sequel, we keep the same notation $∥ \cdot ∥$ to denote the operator norm from $L^{2} (Ω) \to L^{2} (Ω)$ . Using that

(z^{- α} I + T_{h})^{- 1} T_{h} = I - z^{- α} (z^{- α} I + T_{h})^{- 1},

(2.6)

we obtain

∥ (z^{- α} I + T_{h})^{- 1} T_{h} ∥ \leq 1 + M \forall z \in Σ_{θ} .

(2.7)

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

u_{h} (t) = E_{h} (t) u_{0 h} + \int_{0}^{t} {¯ E}_{h} (t - s) f (u_{h} (s)) d s, t > 0.

(2.8)

The operators $E_{h} (t) : L^{2} (Ω) \to L^{2} (Ω)$ and ${¯ E}_{h} (t) : L^{2} (Ω) \to L^{2} (Ω)$ are defined by

E_{h} (t) = \frac{1}{2 π i} \int_{Γ_{θ, δ}} e^{z t} z^{- 1} K_{h} (z) d z and {¯ E}_{h} (t) = \frac{1}{2 π i} \int_{Γ_{θ, δ}} e^{z t} z^{- α} K_{h} (z) d z,

(2.9)

respectively, where $K_{h} (z) := (z^{- α} I + T_{h})^{- 1} T_{h}$ . The contour $Γ_{θ, δ} = {ρ e^{\pm i θ} : ρ \geq δ} \cup {δ e^{i ψ} : | ψ | \leq θ},$ with $θ \in (π / 2, π)$ and $δ > 0$ , is oriented with an increasing imaginary part. Similarly, the solution $u$ of problem (2.1) is given by

u (t) = E (t) u_{0} + \int_{0}^{t} ¯ E (t - s) f (u (s)) d s, t > 0,

(2.10)

where the operators $E$ and $¯ E$ are respectively defined in terms of $K (z) := (z^{- α} I + T)^{- 1} T$ as in (2.9). Standard arguments show that (see for instance [25])

∥ E (t) v ∥ + ∥ E_{h} (t) v ∥ + t^{1 - α} (∥ ¯ E (t) v ∥ + ∥ {¯ E}_{h} (t) v ∥) \leq c ∥ v ∥ \forall v \in L^{2} (Ω) .

(2.11)

Now let $e (t) = u_{h} (t) - u (t)$ denote the error at time $t$ . Define the intermediate solution $v_{h} (t) \in V_{h}$ , $t \geq 0$ , by

T_{h}^{C} \partial_{t}^{α} v_{h} + v_{h} = T_{h} f (u), t > 0, v_{h} (0) = u_{0 h} .

(2.12)

Then, by splitting the error $e = (u_{h} - v_{h}) + (v_{h} - u) =: η + ξ$ , and subtracting (2.12) from (2.1), we find that $ξ$ satisfies

T_{h}^{C} \partial_{t}^{α} ξ (t) + ξ (t) = (T_{h} - T) (f (u) -^{C} \partial_{t}^{α} u) (t), t > 0.

(2.13)

With

ρ (t) := (T_{h} - T) (f (u) -^{C} \partial_{t}^{α} u) (t),

(2.14)

we thus obtain

T_{h}^{C} \partial_{t}^{α} ξ (t) + ξ (t) = ρ (t), t > 0, ξ (0) \in L^{2} (Ω) .

(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 $y$ is a nonnegative function in $L^{1} (0, T_{0})$ which satisfies

y (t) \leq g (t) + β \int_{0}^{t} (t - s)^{- α} y (s) d s for t \in (0, T_{0}],

where $g (t) \geq 0$ , $β \geq 0$ , and $0 < α < 1$ . Then there exists a constant $C_{T_{0}}$ such that

y (t) \leq g (t) + C_{T_{0}} \int_{0}^{t} (t - s)^{- α} g (s) d s for t \in (0, T_{0}] .

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

∥ f (u (t)) ∥ \leq ∥ f (u (t)) - f (0) ∥ + ∥ f (0) ∥ \leq L ∥ u ∥ + ∥ f (0) ∥, t \geq 0,

(2.16)

Lemma 2.1 implies that $∥ u (t) ∥ \leq c (∥ u_{0} ∥ + ∥ f (0) ∥)$ for $t \geq 0$ with $c = c (α, L, T_{0})$ . Now we are ready to prove an error estimate for problem (2.1). Here $ρ$ is given by (2.14) and $~ ρ (t) := \int_{0}^{t} ρ (s) d s$ .

Lemma 2.2.

Let $u$ and $u_{h}$ be the solutions of (2.1) and (2.2), respectively. Assume that $T_{h} (u_{0} - u_{0 h}) = 0$ . Then, for $t > 0$ ,

∥ e (t) ∥ \leq c (G (t) + \int_{0}^{t} (t - s)^{α - 1} G (s) d s),

(2.17)

where

G (t) = t^{- 1} sup s \leq t (∥ ~ ρ (s) ∥ + s ∥ ρ (s) ∥ + s^{2} ∥ ρ_{t} (s) ∥),

(2.18)

and $c$ is independent of $h$ .

Proof.

First we derive a bound for the difference between $u$ and the intermediate solution $v_{h}$ . Since $T_{h} ξ (0) = 0$ , an application of Lemma 3.5 in [16] to (2.1) and (2.12) yields $∥ ξ (t) ∥ \leq c G (t)$ , where $G$ is given by (2.18). In view of the splitting $e = η + ξ$ , it suffices to estimate $∥ η ∥$ . Note that $η$ satisfies

T_{h}^{C} \partial_{t}^{α} η (t) + η (t) = T_{h} (f (u_{h}) - f (u)), t > 0, η (0) = 0.

Hence, by Duhamel’s principle,

η (t) = \int_{0}^{t} {¯ E}_{h} (t - s) (f (u_{h} (s)) - f (u (s))) d s, t > 0.

Using the property of ${¯ E}_{h}$ in (2.11) and condition (1.2), we see that

∥ η (t) ∥ \leq c L \int_{0}^{t} (t - s)^{α - 1} ∥ u_{h} (s) - u (s) ∥ d s, t > 0,

and thus

∥ e (t) ∥ \leq ∥ ξ (t) ∥ + c \int_{0}^{t} (t - s)^{α - 1} ∥ e (s) ∥ d s, 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 $G \in L^{1} (0, T)$ . Recalling that $ρ = (T_{h} - T) L u$ , we have by (2.3), $∥ \partial_{t} ρ (t) ∥ \leq c h^{γ} ∥ L \partial_{t} u (t) ∥$ . Hence, to achieve a $O (h^{γ})$ order of convergence, we need to assume that $\partial_{t} u (t) \in D (L)$ for $t \in (0, T_{0}]$ . 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

S_{h} (z) = (z^{- α} I + T_{h})^{- 1} T_{h} - (z^{- α} I + T)^{- 1} T .

Then $S_{h}$ satisfies the following property.

Lemma 3.1.

There holds

∥ S_{h} (z) v ∥ \leq c h^{γ} | z |^{α (1 - ν)} ∥ L^{ν} v ∥ \forall z \in Σ_{θ}, ν \in [0, 1] .

(3.1)

Proof.

Using the identity (2.6), we verify that

	$S_{h} (z)$	$=$	$z^{- α} (z^{- α} I + T)^{- 1} [(z^{- α} I + T_{h}) - (z^{- α} I + T)] (z^{- α} I + T_{h})^{- 1}$
		$=$	$z^{- α} (z^{- α} I + T)^{- 1} (T_{h} - T) (z^{- α} I + T_{h})^{- 1} .$

Then, by (2.5) and (2.3),

∥ S_{h} (z) v ∥ \leq c ∥ (T_{h} - T) (z^{- α} I + T_{h})^{- 1} v ∥ \leq c h^{γ} | z |^{α} ∥ v ∥ .

This shows (3.1) for $ν = 0$ . For $ν = 1$ , i.e., $v \in D (L)$ , we have $T L v = v$ . Then, by (2.3) and (2.7), we get

∥ S_{h} (z) v ∥ \leq c ∥ (T_{h} - T) (z^{- α} I + T)^{- 1} T L v ∥ \leq c h^{γ} ∥ L v ∥ .

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

We further introduce the following operators: $F_{h} (t) = E_{h} (t) - E (t)$ and ${¯ F}_{h} (t) = {¯ E}_{h} (t) - ¯ E (t)$ . Then, by Lemma 3.1,

∥ {¯ F}_{h} (t) v ∥ \leq c h^{γ} \int_{Γ_{θ, 1 / t}} e^{R e (z) t} | d z | ∥ v ∥ \leq c t^{- 1} h^{γ} ∥ v ∥ .

(3.2)

Similarly, based on (3.1), the following estimate

∥ F_{h} (t) v ∥ \leq c t^{- α (1 - ν)} h^{γ} ∥ L^{ν} v ∥, ν = 0, 1,

(3.3)

holds for $F_{h} (t)$ . Now we are ready to prove a nonsmooth data error estimate. Here and the throughout the paper, $ℓ_{h} (ν) = | ln h |$ if $ν = 0$ and $1$ otherwise.

Theorem 3.1.

Let $u_{0} \in D (L^{ν})$ , $ν \in [0, 1]$ . Let $u$ and $u_{h}$ be the solutions defined by (2.10) and (2.8), respectively, with $u_{0 h} = P_{h} u_{0}$ . Then there is a constant $c = c (r, L, T_{0})$ , where $r \geq ∥ L^{ν} u_{0} ∥ + ∥ f (0) ∥$ , such that

∥ u_{h} (t) - u (t) ∥ \leq c h^{γ} ℓ_{h} (ν) t^{- α (1 - ν)}, t \in (0, T_{0}] .

(3.4)

Proof.

Recall that $e = u_{h} - u$ . From (2.8) and (2.10), we get after rearrangements

e (t) = F_{h} (t) u_{0} + \int_{0}^{t} {¯ E}_{h} (t - s) [f (u_{h} (s)) - f (u (s))] d s + \int_{0}^{t} {¯ F}_{h} (t - s) f (u (s)) d s .

(3.5)

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

I = \int_{0}^{t} {¯ F}_{h} (t - s) (f (u (s)) - f (u (t))) d s and I I = (\int_{0}^{t} {¯ F}_{h} (t - s) d s) f (u (t)) .

For $ν \in (0, 1]$ , we use (3.2), (1.2) and the property $u \in C^{α ν} ([0, T_{0}], L^{2} (Ω))$ to get

∥ I ∥ \leq c h^{γ} \int_{0}^{t} (t - s)^{- 1} (t - s)^{α ν} d s \leq c h^{γ} t^{α ν} .

To estimate $I I$ , we introduce the operator $˜ E (t) = \frac{1}{2 π i} \int_{Γ_{θ, δ}} e^{z t} z^{- 1 - α} S_{h} (z) d z .$ Then ${˜ E}^{'} (t) = ¯ F (t)$ and $∥ ˜ E (t) ∥ \leq c h^{γ}$ for all $t \geq 0$ since $∥ S_{h} (z) ∥ \leq c h^{γ} | z |^{α}$ . Hence, $∥ I I ∥ \leq ∥ f (u (t)) ∥ ∥ ˜ E (t) - ˜ E (0) ∥ \leq c h^{γ}$ as $∥ f (u) ∥$ is bounded, see (2.16). Now, using the properties of ${¯ E}_{h}$ and $F_{h}$ in (2.11) and (3.3), respectively, we obtain

∥ e (t) ∥ \leq c t^{- α (1 - ν)} h^{γ} ∥ L^{ν} u_{0} ∥ + c \int_{0}^{t} (t - s)^{α - 1} ∥ e (s) ∥ d s + c h^{γ} + c h^{γ} t^{α ν} .

(3.6)

The desired estimate follows now by applying Lemma 2.1. To establish the estimate for $ν = 0$ , 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 $V_{h} \subset H_{0}^{1} (Ω)$ consists of piecewise linear functions on a shape-regular triangulation with a mesh parameter $h$ . Let $T_{h} : L^{2} (Ω) \to V_{h}$ be the solution operator of the discrete problem:

a (T_{h} f, v) = (f, v) \forall v \in V_{h},

(3.7)

where $a (\cdot, \cdot)$ is the bilinear form associated with the elliptic operator $L$ . Then, $T_{h}$ is selfadjoint, positive semidefinite on $L^{2} (Ω)$ and positive definite on $V_{h}$ , see [3], and satisfies (2.3) with $γ = 2$ . Thus, by Theorem 3.1,

∥ u_{h} (t) - u (t) ∥ \leq c h^{2} ℓ_{h} (ν) t^{- α (1 - ν)}, t > 0,

(3.8)

for $ν \in [0, 1]$ . This improves the following estimate

∥ u_{h} (t) - u (t) ∥ \leq c h^{2} (t^{- α (1 - ν)} + max (0, ln (t^{α (1 - ν)} / h^{2}))), t > 0,

(3.9)

established in [1]. We notice that the logarithmic factor is also present in the parabolic case when $ν = 0$ , 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 (\partial_{t}) u + L u = f (u) in Ω \times (0, T_{0}], u (0) = u_{0} in Ω, u = 0 on \partial Ω \times (0, T_{0}],

(3.10)

where the multi-term differential operator $P (\partial_{t})$ is defined by $P (\partial_{t}) = \partial_{t}^{α} + \sum_{i = 1}^{m} b_{i} \partial_{t}^{α_{i}}$ with $0 < α_{m} \leq \dots \leq α_{1} \leq α < 1$ being the orders of the fractional Caputo derivatives, and $b_{i} > 0$ , $i = 1, \dots, m$ . 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 $F_{h}$ , ${¯ E}_{h}$ and ${¯ F}_{h}$ . 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 $α = 1$ .

Remark 3.2.

For smooth initial data $u_{0} \in D (L)$ , the estimate (3.4) still holds for the choice $u_{0 h} = R_{h} u_{0}$ . Indeed, we have

E_{h} (t) R_{h} u_{0} - E (t) u_{0} = E_{h} (t) (R_{h} u_{0} - u_{0}) + (E_{h} (t) u_{0} - E (t) u_{0}) .

By the stability of the operator $E_{h} (t)$ ,

∥ E_{h} (t) (R_{h} u_{0} - u_{0}) ∥ \leq c ∥ R_{h} u_{0} - u_{0} ∥ \leq c h^{γ} ∥ L u_{0} ∥ .

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 $[0, T_{0}]$ into $N$ equal subintervals with a time step size $τ = T_{0} / N$ , and let $t_{j} = j τ$ . The convolution quadrature [23] refers to an approximation of any function of the form $k * φ$ as

(k * φ) (t_{n}) := \int_{0}^{t_{n}} k (t_{n} - s) φ (s) d s \approx n \sum j = 0 β_{n - j} (τ) φ (t_{j}),

where the weights $β_{j} = β_{j} (τ)$ are computed from the Laplace transform $K (z)$ of $k$ rather than the kernel $k (t)$ . With $\partial_{t}$ being time differentiation, define $K (\partial_{t})$ as the operator of (distributional) convolution with the kernel $k$ : $K (\partial_{t}) φ = k * φ$ for a function $φ (t)$ with suitable smoothness. Then a convolution quadrature will approximate $K (\partial_{t}) φ$ by a discrete convolution $K (\partial_{τ}) φ$ at $t = t_{n}$ as $K (\partial_{τ}) φ (t_{n}) = \sum_{j = 0}^{n} β_{n - j} (τ) φ (t_{j}),$ where the quadrature weights ${β_{j} (τ)}_{j = 0}^{\infty}$ are determined by the generating power series $\sum_{j = 0}^{\infty} β_{j} (τ) ξ^{j} = K (δ (ξ) / τ)$ 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, $δ (ξ) = 1 - ξ$ .

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

K_{2} (\partial_{τ}) K_{1} (\partial_{τ}) = K_{2} K_{1} (\partial_{τ}) and K_{2} (\partial_{τ}) (k_{1} * φ) = (K_{2} (\partial_{τ}) k_{1}) * φ .

(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 $G (z)$ be analytic in the sector $Σ_{θ}$ and such that

∥ G (z) ∥ \leq c | z |^{- μ} \forall z \in Σ_{θ},

for some real $μ$ and $c$ . Then, for $φ (t) = c t^{σ - 1}$ , the convolution quadrature based on the backward Euler method satisfies

(4.2)

Upon using the relation between the Riemann-Liouville derivative denoted by $\partial_{t}^{α}$ and the Caputo derivative $^{C} \partial_{t}^{α}$ , the semidiscrete scheme (2.2) can be rewritten as

T_{h} \partial_{t}^{α} (u_{h} - u_{0 h}) + u_{h} = T_{h} f (u_{h}), t > 0, u_{h} (0) = u_{0 h} .

(4.3)

Thus, the proposed backward Euler CQ scheme is to seek $U_{h}^{n} \in V_{h}$ , $n \geq 0$ , such that

T_{h} \partial_{τ}^{α} (U_{h}^{n} - U_{h}^{0}) + U_{h}^{n} = T_{h} f (U_{h}^{n}), n \geq 1, U_{h}^{0} = u_{0 h} .

(4.4)

4.1 The linear case

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

T^{C} \partial_{t}^{α} u + u = T f (t), t > 0, u (0) = u_{0},

(4.5)

with a semidiscrete solution $u_{h} (t) \in V_{h}$ satisfying

T_{h}^{C} \partial_{t}^{α} u_{h} + u_{h} = T_{h} f (t), t > 0, u_{h} (0) = P_{h} u_{0} .

(4.6)

The fully discrete solution $U_{h}^{n} \in V_{h}$ is defined by

T_{h} \partial_{τ}^{α} (U_{h}^{n} - U_{h}^{0}) + U_{h}^{n} = T_{h} f (t_{n}), n \geq 1, U_{h}^{0} = P_{h} u_{0} .

(4.7)

Then we establish the following result.

Theorem 4.1.

Let $u_{0} \in D (L^{ν})$ , $ν \in [0, 1]$ . Let $u$ and $U_{h}^{n}$ be the solutions of $(???)$ and $(???)$ , respectively, with $f = 0$ . Then, for $t_{n} > 0$ ,

∥ U_{h}^{n} - u (t_{n}) ∥ \leq c (τ t_{n}^{α ν - 1} + h^{γ} t_{n}^{- α (1 - ν)}) ∥ L^{ν} u_{0} ∥ .

(4.8)

Proof.

We first notice that $∥ u_{h} (t_{n}) - u (t_{n}) ∥ \leq c t_{n}^{- α (1 - ν)} h^{γ} ∥ L^{ν} u_{0} ∥$ , which follows from Theorem 3.1 in the case $f = 0$ . In order to estimate $∥ U_{h}^{n} - u_{h} (t_{n}) ∥$ , apply $\partial_{t}^{- α}$ and $\partial_{τ}^{- α}$ to (4.6) and (4.7), respectively, use the associativity of convolution and the property $T_{h} P_{h} = T_{h}$ , to deduce that

U_{h}^{n} - u_{h} (t_{n}) = - (G (\partial_{τ}) - G (\partial_{t})) u_{0},

(4.9)

where $G (z) = (z^{- α} I + T_{h})^{- 1} T_{h}$ . We recall that, by (2.5), $∥ G (z) ∥ \leq c \forall z \in Σ_{θ} .$ Then, Lemma 4.1 (with $μ = 0$ , $σ = 1$ ) and the $L^{2}$ -stability of $P_{h}$ yield

∥ U_{h}^{n} - u_{h} (t_{n}) ∥ \leq c τ t_{n}^{- 1} ∥ u_{0} ∥ .

(4.10)

For $u_{0} \in D (L)$ , consider first the choice $u_{h} (0) = R_{h} u_{0}$ . Recalling that $R_{h} = T_{h} L$ , we use the identity $G (z) = I - z^{- α} (z^{- α} I + T_{h})$

Galerkin Type Methods for Semilinear Time-Fractional Diffusion Problems

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

Numerical Analysis of Backward Subdiffusion Problems

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

Complete monotonicity-preserving numerical methods for time fractional ODEs

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

Incomplete Iterative Solution of the Subdiffusion Problem

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

1 Introduction

2 General setting and preliminary estimates

Lemma 2.1.

Lemma 2.2.

Proof.

3 Error estimates without regularity assumptions

Lemma 3.1.

Proof.

Theorem 3.1.

Proof.

Remark 3.1.

Remark 3.2.

4 Fully discrete schemes

Lemma 4.1.

4.1 The linear case

Theorem 4.1.

Proof.

Galerkin Type Methods for Semilinear Time-Fractional Diffusion Problems

Related Research

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

Numerical Analysis of Backward Subdiffusion Problems

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

Complete monotonicity-preserving numerical methods for time fractional ODEs

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

Incomplete Iterative Solution of the Subdiffusion Problem

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

1 Introduction

2 General setting and preliminary estimates

Lemma 2.1.

Lemma 2.2.

Proof.

3 Error estimates without regularity assumptions

Lemma 3.1.

Proof.

Theorem 3.1.

Proof.

Remark 3.1.

Remark 3.2.

4 Fully discrete schemes

Lemma 4.1.

4.1 The linear case

Theorem 4.1.

Proof.