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
(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
(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
(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 domains . Note that
form a Hilbert scale of interpolation spaces and
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]:
(1.4) |
(1.5) |
(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])
(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
(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
(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
(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
(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
(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
(2.6) |
we obtain
(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
(2.8) |
The operators and are defined by
(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
(2.10) |
where the operators and are respectively defined in terms of as in (2.9). Standard arguments show that (see for instance [25])
(2.11) |
Now let denote the error at time . Define the intermediate solution , , by
(2.12) |
Then, by splitting the error , and subtracting (2.12) from (2.1), we find that satisfies
(2.13) |
With
(2.14) |
we thus obtain
(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
where , , and . Then there exists a constant such that
Note that, by using (2.10), (2.11) and the inequality
(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.
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
Hence, by Duhamel’s principle,
Using the property of in (2.11) and condition (1.2), we see that
and thus
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
Then satisfies the following property.
Lemma 3.1.
There holds
(3.1) |
Proof.
We further introduce the following operators: and . Then, by Lemma 3.1,
(3.2) |
Similarly, based on (3.1), the following estimate
(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.
Proof.
Recall that . From (2.8) and (2.10), we get after rearrangements
(3.5) |
The last term in (3.5) can be written as where
For , we use (3.2), (1.2) and the property to get
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
(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:
(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,
(3.8) |
for . This improves the following estimate
(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:
(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.
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
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
(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
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
(4.3) |
Thus, the proposed backward Euler CQ scheme is to seek , , such that
(4.4) |
4.1 The linear case
We begin by investigating the time discretization of the inhomogeneous linear problem
(4.5) |
with a semidiscrete solution satisfying
(4.6) |
The fully discrete solution is defined by
(4.7) |
Then we establish the following result.
Theorem 4.1.
Let , . Let and be the solutions of and , respectively, with . Then, for ,
(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
(4.9) |
where . We recall that, by (2.5), Then, Lemma 4.1 (with , ) and the -stability of yield
(4.10) |
For , consider first the choice . Recalling that , we use the identity