1. Introduction
We analyze a new method for numerically solving a spectral fractional diffusion equation in a curvilinear polygonal domain
. There are no additional constraints on the geometry of
. The equation can be written in the form
, where
is a positive definite self-adjoint second-order elliptic operator in
with homogeneous Dirichlet data. In recent years, several groups of methods have been proposed to approximate the inverse of
in the case of sub diffusion, i.e., for
. Despite the different approaches in obtaining them, they can all be written as a rational approximation. In the original BURA (Best Uniform Rational Approximation) method, we used the approximation of
,
. In the following years, this approach was applied to more general equations that involve
. By definition, when applicable, BURA methods have the best accuracy. The same applies to their computational complexity as well. Let us consider the second-order elliptic equation with homogeneous Dirichlet data
Here,
is a bounded curvilinear polygonal domain in
, and we assume that
for
. The fractional powers of the elliptic operator associated with the problem (
1) are defined in terms of the weak form of the differential equation, namely,
is the unique function in
satisfying
where
For
, (
2) defines a solution operator
. Following (
1), we introduce the unbounded operator
with domain
such that
for
where
with
.
is well defined as
is injective.
Now, we consider the spectral fractional diffusion problem
where
Here,
is the fractional power of the diffusion operator,
are the normalized eigenfunctions of
, satisfying the equalities
, also assuming that the positive eigenvalues
are in monotonically increasing ordering, i.e.,
if
.
We apply the finite element method (FEM) to the numerical solution of (
3), thus obtaining the linear system
Here,
is a sparse symmetric positive definite (SPD) matrix,
where
is the diagonal matrix of the positive eigenvalues
of
, and
is the matrix of the corresponding eigenvectors
. Analogous to the continuous case, here we assume that the eigenvectors of
are normalized to satisfy the equalities
, and
if
.
In the case of FEM discretization, the matrix
is symmetric with respect to the energy scalar product associated with the mass matrix. The details are discussed in
Section 2.
The Formulas (
4) and (
6) are clear and simple to interpret. But, are they applicable in computational practice? In this regard, the first question is if the spectrum of
or
is known, or if it is easily computable. In the general case (regarding the geometry of
in particular), this is not the case. And, in the event that the spectrum is available, the problem of implementation remains. In the discrete case, this concerns the computational complexity to perform matrix-vector multiplications with
and
, which, excluding some very special cases, is of order
. Therefore, the relations (
4) and (
6) can be mainly used in theoretical analyses.
There are different definitions of fractional Laplacian operators, and more generally of fractional diffusion operators. They correspond to nonequivalent boundary value problems. Fractional diffusion operators are non-local. As a result, the matrices emerging after their finite element discretization are dense. One of the reasons for the growing interest in spectral fractional diffusion is the recognized opportunities for development of efficient numerical methods for multidimensional problems. Advances on this topic have been strongly influenced by the work of Caffarelli and Silvestre [
1] (see also [
2]) on an extension problem associated with the fractional Laplacian. Almost ten years later, this result has been used to develop a numerical method [
3] where the non-locality is avoided at the cost of the increased dimensions of the computational domain.
Transformations of the boundary value problem with fractional Laplacian in
to
-dimensional local (of non-fractional order) equations have been used in a number of approaches proposed over the past decade. From this view point, and without pretending for completeness, we will cite the following classification used, for example, in [
4]:
- AP1.
Methods based on the elliptic extension in
, see, e.g., [
3,
5,
6];
- AP2.
Methods based on the pseudo-parabolic extension in
, see, e.g., [
7,
8,
9];
- AP3.
Methods based on the Dunford–Taylor integral representation of the solution, see, e.g., [
10,
11];
- AP4.
Methods based on the best uniform rational approximation (BURA), see, e.g., [
12,
13].
Despite the substantial difference of the approaches AP1–AP4, the implementation of all of them can be interpreted as some rational approximation of
. This common property was first systematically discussed in [
14]; see also the survey paper [
15]. Over the next few years, the work on the further development of efficient numerical methods for space-fractional differential equations is actively continuing [
4,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28].
In the recent work [
20], the convergence of
FEM for spectral fractional diffusion in polygons is studied. The developed method is based on the elliptic extension in
, i.e., on the approach AP1. Following the spectral decomposition proposed in [
5], the extended elliptic operator is diagonalized with respect to the auxiliary variable in
. Then,
FEM is applied to the obtained series of auxiliary singularly perturbed diffusion–reaction problems in
. Using the error analysis published in [
29], exponential convergence is proven.
Our research is motivated by [
20]. However, the construction of the proposed method is significantly different. The goal is to improve the computational complexity as well as reduce some numerical stability limitations associated with the operator diagonalization used in [
20]. While ref. [
20] exploits the elliptic extension of Caffarelli and Silvestre [
1], this new method is based on best uniform rational approximation of the solution. The similarities/differences and advantages/disadvantages of these two approaches to numerically solving fractional diffusion problems are discussed in detail in the survey paper [
12]. The main contributions of this paper are as follows. The developed and analyzed semi-discrete BURA method is substantially new. The abbreviation BURA-SD is used to distinguish the new method from the previously known fully discrete BURA method. The semi-discrete scheme allows for the better adaptation of the construction and analysis of the rational approximation to the regularity of the solution in the case of a polygonal domain. To maximize the quality of the overall computational complexity, we integrate the optimal exponential convergence rate of BURA-SD and the exponential convergence rate of
FEM.
The rest of the paper is organized as follows. Basic results about
FEM discretization in curvilinear polygonal domains are presented in
Section 2.
Section 3 discusses how to combine BURA-SD and
FEM discretization. The next two sections contain the main results of the paper. Exponential error estimates of the semi-discrete and fully discrete BURA-SD methods are obtained in
Section 4. The computational complexity is then analyzed by balancing the different types of errors in the next section. At the end, brief concluding remarks and challenging topics for future research are given.
2. FEM Discretization
The finite element approximation is defined in terms of a conforming finite dimensional space
of piece-wise polynomial functions over a triangulation
of the polygonal domain
. Here, we consider isoparametric Lagrangian finite elements (FE) of degree
. Following the assumptions of [
20], we allow both triangular and quadrilateral elements
, which in the general case for
can be curvilinear. In principle, shape regularity of the triangulation is not required. In fact, as we shall discuss later, anisotropic mesh refinement towards
is natural to be applied to resolve singularities of the boundary layers or around corner points of the boundary. Such singularities are generically present in solutions of fractional diffusion problems in space.
The discrete operator
is defined as the inverse of
with
, where
is the unique solution to
Then, the finite element numerical solution of (
3), i.e., the approximation
of
u, is defined by the equation
Here,
stands for the
projection onto
. In this case,
N denotes the dimension of the FEM space
and equals the number of (interior with respect to
) degrees of freedom. The operator
in the finite element case is a map of
into
so that
, where
is the unique solution to
Let us denote by
the standard nodal basis of
. With respect to this basis, the operator
is represented by the matrix
In accordance with the terminology adopted in the finite element method,
and
are the mass and stiffness matrices, respectively. They are both symmetric and positive definite (SPD).
Obviously, if
and
are the coefficient vectors corresponding to
, then
. Now, for the coefficient vector
corresponding to
, we obtain
, where
is the vector with entries
Thus, using vector notation so that
is the coefficient vector representing the solution
through the nodal basis, we can write the finite element approximation of (
3) in the form of a system
Therefore, for the finite element approximation of the subdiffusion problem (
7) we obtain the equation
We note that the matrix
is SPD with respect to the Euclidean space defined by the energy inner product associated with the mass matrix, that is
This follows directly from the definition (
10), which implies the equality
In the context of our study, the FEM problem (
8) is by default large-scale, i.e.,
. When analyzing the new BURA-SD algorithm, we will assume that systems of the type
,
and
can be solved numerically in an efficient way by a method that requires
arithmetic operations, with a constant possibly depending on
p. This can be achieved by using fast iterative solvers for the aforementioned auxiliary sparse linear systems, based for example on multigrid, multilevel or domain decomposition preconditioners, see, e.g., [
30]. In the present paper, our goal is to construct a solution method for (
3) whose total computational complexity is as close to
as possible.
In the case of quasi-uniform partition
and
h FEM approximation, the condition number of the stiffness matrix
increases as
when
, i.e.,
. Here, the focus is on the
FEM. It is known that in this case, shape functions based on interpolating Lagrange polynomials at Gauss–Lobato points play an important role. Now, let us denote by
the standard nodal basis stiffness matrix and by
the stiffens matrix obtained if the corresponding spectral basis functions are used in the interior of the finite elements
. If the mesh
is quasi-uniform, then the following estimates hold (see [
29]):
As mentioned above, it is natural to apply local mesh refinement towards
to resolve singularities of the boundary layers or around corner points of the boundary. The case of an L-shaped domain with a mesh refined towards the re-entrant corner is discussed in [
31], see also [
32]. From the presented numerical test, one can see in particular the exponential growth of
when the number of layers increases. The behavior of
is significantly better. Under rather more general conditions, it was shown in [
29] that, for meshes that are refined geometrically toward singularities,
is independent of the number of refinement levels.
The condition number of the mass matrix is also of interest when the BURA
CD solver is implemented. In the case of quasi-uniform partition
and
h FEM approximation,
. Now, unlike (
13), if
is quasi-uniform and the
FEM is applied,
These estimates do not depend on
h, which is a general property of the condition number of the mass matrix. Since the element mass matrix is positive definite, estimates of this type are proved element-wise. And again, under conditions similar to those in [
29], for meshes that are geometrically refined in the singularity region,
does not depend on the number of levels of refinement.
Estimates (
13) and (
14) are important for a better understanding of the computational complexity analysis in
Section 5.
Finally, although the improved BURA method presented in [
12] is robust with respect to the condition number of the discrete elliptic operator, the estimate of
is of practical importance in computing the best uniform rational approximation. So, for example, in the numerical tests, we use the BRASIL [
33] software to calculate BURA, where an upper bound of
is needed as an input parameter.
3. BURA-SD
Let us denote by
the best uniform rational approximation (BURA) of degree
k of
in
, belonging to the diagonal of the Walsh table. This means that we consider the rational functions
, where
and
are polynomials of equal degree
k. By definition,
is the error minimizer for which
The following additive and multiplicative representations of
can be used in the implementation of the BURA methods and algorithms:
As was shown in [
13], all
k zeros
and poles
are real and negative. Even more precisely, they are interlacing, that is, with appropriate numbering
The following asymptotically sharp exponential estimate of
(see for more details [
13] and references there in) holds:
In accordance with the introduced notations, the BURA approximation of
for
is defined as
where
is the smallest eigenvalue of
. In the particular case of FEM approximation, see (
10),
is the smallest positive eigenvalue of the generalized eigenvalue problem
Here, we will use the variant of the BURA method as proposed and studied in [
12], see also the survey paper [
15]. Let us denote by
the BURA approximation of the solution
of
FEM system (
11) defined as follows
Following (
16), we distinguish additive and multiplicative variants of the BURA method. So far, the first one is more often used and analyzed. In this paper, we will look at both. The additive variant relies on the partial fraction decomposition of
, while the multiplicative one deals with the factorization of
. Thus, transforming the expressions (
16), we obtain the equalities
Here,
and
. In this way, the variants of BURA method based on (
19) read as
Thus, the additive and multiplicative variants of the BURA method for numerical solution of (
3) require solving
k auxiliary linear systems with SPD matrices, which are positive diagonal shifts of
. We should note here that the matrix
is dense and we do not want to compute it explicitly and then solve the systems that appear in (
20). To overcome this problem we rewrite the additive and multiplicative representations (
20) in the form
Thus, the determining part of the computational complexity the BURA method is due to solving
k linear systems with sparse SPD matrices, which are positive shifts of the stiffness matrix
with the mass matrix
.
The discussed BURA method can be directly applied to the
FEM discretization of (
3) as follows:
- (1)
Consider the fractional diffusion problem (
3),
,
.
- (2)
Choose the FEM space and compute the stiffness matrix corresponding to the elliptic operator and the mass matrix .
- (3)
Choose the degree of rational approximation k and find the best uniform rational approximation of , .
- (4)
Applying (
18), compute the BURA approximation of the solution
of the system (
11),
.
- (5)
Find the BURA approximation of as the FEM function represented by the nodal basis vector .
We now introduce the BURA-SD method. In the new method, the best uniform rational approximation is applied to the continuous fractional diffusion equation. As we will see, the semi-discrete BURA approximation leads to a composite method whose implementation and analysis better suits the problem at hand.
- (1)
Consider the fractional diffusion problem (
3),
,
.
- (2)
Choose the degree of rational approximation k and find the best uniform rational approximation of , .
- (3)
Find the semi-discrete BURA approximation
of
u by the formula
- (4)
Choose the
FEM space
and compute the stiffness matrix
corresponding to the elliptic operator
and the mass matrix
. Apply the
FEM to (
22) and obtain the approximation
of
, which satisfies the equality
- (5)
Find the BURA approximation of as the FEM function represented by the nodal basis vector .
Remark 1. The BURA and BURA-SD methods are not equivalent, but the respective approximations asymptotically coincide. More precisely, it is easy to see that Let us recall that the first eigenvalue of the operator associated with the elliptic boundary value problem (1) is uniformly separated from 0. Thus, and accordingly, for a sufficiently fine FEM discretization, . We will also note that in the implementation of both methods, BURA and BURA-SD, some approximation from below of the first eigenvalue is used. In numerically solving spectral fractional diffusion problems, it is generally assumed that the power
k required to obtain a given accuracy of
is a measure of computational efficiency. This principle will also be applied in the analysis in
Section 5. However, it is worth noting that some recent publications raise questions and suggest ideas for further improvements in this direction as well. Thus, for example, it was shown in [
27] that for larger
k, the so-called additive reduced sum BURA-AR method and the multiplicative reduced product BURA-MR method can be applied to improve the computational complexity. Although the approach is different, the results based on the reduced conjugate gradient method proposed in [
34] lead to quite similar conclusions.
5. Computational Complexity
In this section, we investigate sufficient conditions for balancing the terms on the right-hand side of the error estimate in Theorem 1. This is the condition for ensuring the high performance of the proposed method. In (
38), the square root of
k exponential convergence of the BURA-SD interacts with the
FEM exponential accuracy with respect to the mesh parameters
, which are in turn controlled by the scale resolution condition.
An inherent property of the BURA methods is the exponential clustering of the poles and roots of the corresponding best uniform rational approximation [
36], see also [
13,
37]. Here, we will use the inequality
with positive constants
and
that can only depend on the fractional power
.
The exponential decrease of
with increasing
i for
is shown in
Figure 2, where the influence of
is also illustrated. In a local neighborhood of
k, we note a zone of faster decrease. However, this feature is not essential for the estimation of
, taking into account the global monotonicity of
. The BRASIL software [
33] is used to calculate
, see also [
38].
Lemma 4. LetThen, the scale resolution condition (37) is satisfied. Proof. We write (
40) in the form
Then
and therefore
Thus, from (
41) and (
39) we obtain that (
37) is satisfied. □
Corollary 1. The scale resolution condition requires that the number of levels L in the geometrically refined mesh be large enough to provide the necessary FEM accuracy to the -perturbed reaction-diffusion problems that arise in the BURA-SD method. It follows from Lemma 4 that there exists a constant that can only depend on α, so the scale resolution condition is satisfied if We will now derive a sufficient condition for the degree
k of BURA-SD to guarantee balancing of the terms on the right-hand side of the error estimate (
38). First, the BURA-SD and
h FEM errors are balanced if
Therefore,
and consequently
Analogously, the sufficient condition for balancing the BURA-SD and
p FEM errors has the form
Then,
and consequently
Using (
43) and (
44) we obtain the following lemma.
Lemma 5. Let the scale resolution condition (37) be satisfied and letThen,Similar to Lemma 3, the constants depend only on the data and on the macrotriangulation . We conclude our analysis with the following theorem.
Theorem 2. Let be a curvilinear polygon with mesh constructed according to assumptions (A1–A5). Then, there exist a mesh refinement parameter , BURA-SD degree k, and constants , , , (depending only on the data , , and on the macro-triangulation ), such that the following error estimate holdsand the computational complexity is bounded by Proof. The error estimate (
45) follows directly from Lemma 5. Now, to prove the computational complexity estimate (
46) we analyze the implementation of the BURA-SD method. It involves solving
k linear systems with the sparse SPD matrices
,
. We will assume that a suitable fast preconditioned iterative solver is used for this purpose. For example, such a multigrid preconditioner for the
p-hierarchical basis finite element method is proposed in [
30]. The number of multigrid iterations
is uniformly bounded regardless of
p. Then, the following estimate of the computational complexity
of the BURA-SD method holds
where
represents the number of nonzero off-diagonal elements of the matrices. Thus, from (
32) we obtain
and finally applying (
42) we obtain
□
6. Concluding Remarks
In this paper, we have developed a new method for the numerical solution of sub-diffusion problems in curvilinear polygons . This means that there are practically no restrictions on the geometry of the bounded computational domain. The spectral definition of the elliptic operator , , is used. The proposed methodology allows the construction of a method that fully utilizes the strengths of the FEM and new BURA-SD methods. Only less than ten years ago, we did not have the theoretical basis for creating computationally efficient methods and algorithms for large-scale problems of this type.
The error estimate (
45) fully reproduces the exponential convergence rate
of the
FEM, see (
31), which characterizes the case of a standard (local) elliptic boundary value problem. In addition, the developed algorithm providing this accuracy is computationally highly efficient. In this context, we note that taking into account (
32), we can write (
46) in the compact form
. We recall that in constructing the
FEM patches, we denote by
L the number of geometric refinement layers to an edge where
n is the number of geometric refinement layers to a corner that is a vertex of the same edge.
Our research is influenced and motivated by the recently published results in [
20]. In short, the following steps apply there. The auxiliary elliptic extension equation in
that was originally proposed in [
1] is considered. It is first truncated and diagonalized in the expanded variable. The fractional diffusion problem thus reduces to
decoupled singularly perturbed diffusion-reaction equations in
. The
FEM discretization and error analysis from [
35] are then utilized. By the arguments of [
14], we can interpret the method of [
20] as a rational approximation of
. In the context of the present paper, the integer
corresponds to the BURA-SD degree
k.
In this work, we elaborated a novel approach to the BURA (here called BURA-SD) approximation of that is adapted to the specifics in the analysis of the FEM used. One may wonder whether some other rational approximation can be integrated into this framework instead of BURA-SD. The answer is yes. Alternatively, both pseudo-parabolic extension (AP2 approach) and quadrature integral representation formulas (AP3 approach) can be applied. It is also worth noting that in each of the AP2-AP3 cases, the problem of computational instability that appears at the disorganization step of A1 approach for larger is overcome. Compared to all others, the essential advantage of the BURA-based method is that it provides, by definition, the best uniform rational approximation. When applicable, this leads directly to a better computational complexity as well.
The application of a newly developed numerical method for a class of elliptic problems to the corresponding parabolic equations is usually expected to be a direct consequence. Although this is true for standard (local) partial differential equations, this does not hold in the case of space fractional diffusion problems. The reason for this is that multiplying a vector by the matrix
turns out to be more complicated than solving a system with the same matrix. A systematic analysis of BURA and BURA-based approximations of the fractional powers of sparse SPD matrices was recently presented in [
4]. It follows from the results in [
4] that the application of the BURA-SD method proposed here for the case of the
FEM space discretization to parabolic equations will require a thorough analysis of the computational complexity, which will also include the influence of the condition number
.
The degree of rational approximation
k determines the number of auxiliary sparse linear systems to solve in the implementation of the numerical methods for fractional diffusion equations, which we classified here into four groups according to approaches AP1–AP4. This is the argument for using
k as a universal measure when estimating computational complexity via
. More recently, a reduced conjugate gradient basis method for fractional diffusion was proposed in [
34]. It is shown there that a smaller number of linear systems can be used without the loss of accuracy of the rational approximation. This is also the aim of the earlier work [
26], where the so-called reduced multiplicative (BURA-MR) and additive (BURA-AR) methods were presented. Further theoretical analysis is needed to better understand the possible similarities of the apparently quite different approaches explored in [
26,
34]. In any case, such kinds of results can be used to improve the algorithm and thus to refine the estimate of the computational complexity (
46).
The results in [
20], and now in the present paper, raise a wide range of questions for future research. Here, we also note the recent work [
39], where an a posteriori error estimator for the spectral fractional power of the Laplacian is developed. In this context, we will note such challenging examples as the generalization of currently available results to fractional diffusion problems in three-dimensional polygons, as well as to equations with heterogeneous and anisotropic coefficients. The development of a BURA-based multilevel adaptive scheme is another attractive topic for future research.
In conclusion, the discussions in this paper have been focused on mesh methods and in particular on the FEM discretization case. In perspective, however, we are analyzing alternative possibilities to further develop our research in the spirit of meshless methods.