Next Article in Journal
Building an Analog Circuit Synapse for Deep Learning Neuromorphic Processing
Previous Article in Journal
Analyzing Treatment Effect by Integrating Existing Propensity Score and Outcome Regressions with Heterogeneous Covariate Sets
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Exponential Convergence and Computational Efficiency of BURA-SD Method for Fractional Diffusion Equations in Polygons

by
Svetozar Margenov
Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
Mathematics 2024, 12(14), 2266; https://doi.org/10.3390/math12142266
Submission received: 1 July 2024 / Revised: 16 July 2024 / Accepted: 17 July 2024 / Published: 19 July 2024
(This article belongs to the Topic Numerical Methods for Partial Differential Equations)

Abstract

:
In this paper, we develop a new Best Uniform Rational Approximation-Semi-Discrete (BURA-SD) method taking into account the singularities of the solution of fractional diffusion problems in polygonal domains. The complementary capabilities of the exponential convergence rate of BURA-SD and the h p FEM are explored with the aim of maximizing the overall performance. A challenge here is the emerging singularly perturbed diffusion–reaction equations. The main contributions of this paper include asymptotically accurate error estimates, ending with sufficient conditions to balance errors of different origins, thereby guaranteeing the high computational efficiency of the method.

1. Introduction

We analyze a new method for numerically solving a spectral fractional diffusion equation in a curvilinear polygonal domain Ω R 2 . There are no additional constraints on the geometry of Ω . The equation can be written in the form A α u ( x ) = f ( x ) , where A 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 A α in the case of sub diffusion, i.e., for α ( 0 , 1 ) . 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 z α , z [ 0 , 1 ] . In the following years, this approach was applied to more general equations that involve A α . 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
· ( a ( x ) v ( x ) ) = g ( x ) , for x Ω , v ( x ) = 0 , for x Ω .
Here, Ω is a bounded curvilinear polygonal domain in R 2 , and we assume that 0 < a 0 a ( x ) for x Ω . 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, v ( x ) is the unique function in V = H 0 1 ( Ω ) = { w ( x ) H 1 ( Ω ) : w ( x ) = 0 for x Ω } satisfying
A ( v , w ) = ( g , w ) , w V ,
where
A ( v , w ) : = Ω a ( x ) v ( x ) · w ( x ) d x and ( v , w ) : = Ω v ( x ) w ( x ) d x .
For g X : = L 2 ( Ω ) , (2) defines a solution operator Y g : = v . Following (1), we introduce the unbounded operator A with domain
D ( A ) = { Y g : g X }
such that A v = g for u D ( A ) where g X with Y g = u . A is well defined as Y is injective.
Now, we consider the spectral fractional diffusion problem
A α u = f , with   a   solution u = A α f , α ( 0 , 1 ) ,
where
A α u = j = 1 λ j α ( u , ϕ j ) ϕ j , and   therefore u = j = 1 λ j α ( f , ϕ j ) ϕ j .
Here, α is the fractional power of the diffusion operator, { ϕ j } j = 1 are the normalized eigenfunctions of A , satisfying the equalities ( A ϕ i , ϕ j ) = ( ϕ i , ϕ j ) = δ i , j , also assuming that the positive eigenvalues { λ j } j = 1 are in monotonically increasing ordering, i.e., 0 < λ i λ j if i < j .
We apply the finite element method (FEM) to the numerical solution of (3), thus obtaining the linear system
A α u h = f h .
Here, A R N × N is a sparse symmetric positive definite (SPD) matrix,
A α = W t D α W , and   consequently A α = W t D α W ,
where D is the diagonal matrix of the positive eigenvalues { λ i , h } i = 1 N of A , and W is the matrix of the corresponding eigenvectors { ϕ i , h } i = 1 N . Analogous to the continuous case, here we assume that the eigenvectors of A are normalized to satisfy the equalities ( A ϕ i , h , ϕ j , h ) = ( ϕ i , h , ϕ j . h ) = δ i , j , and 0 < λ i , h λ j , h if i < j .
In the case of FEM discretization, the matrix A 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 A or A 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 W and W t , which, excluding some very special cases, is of order O ( N 2 ) . 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 Ω R d to ( d + 1 ) -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 Ω × ( 0 , ) , see, e.g., [3,5,6];
AP2.
Methods based on the pseudo-parabolic extension in Ω × ( 0 , 1 ) , 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 A α . 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 h p FEM for spectral fractional diffusion in polygons is studied. The developed method is based on the elliptic extension in Ω × ( 0 , ) , 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 ( 0 , ) . Then, h p 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 h p FEM.
The rest of the paper is organized as follows. Basic results about h p FEM discretization in curvilinear polygonal domains are presented in Section 2. Section 3 discusses how to combine BURA-SD and h p 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. hp FEM Discretization

The finite element approximation is defined in terms of a conforming finite dimensional space V h , p V of piece-wise polynomial functions over a triangulation T h of the polygonal domain Ω . Here, we consider isoparametric Lagrangian finite elements (FE) of degree p 1 . Following the assumptions of [20], we allow both triangular and quadrilateral elements K T h , which in the general case for p > 1 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 A h is defined as the inverse of Y h : V h , p V h , p with Y h g h : = v h , where v h V h , p is the unique solution to
A ( v h , w h ) = ( g h , w h ) , for   all w h V h , p .
Then, the finite element numerical solution of (3), i.e., the approximation u h V h , p of u, is defined by the equation
A h α u h = π h f , or   equivalently u h = A h α π h f : = A h α f h .
Here, π h stands for the L 2 ( Ω ) projection onto V h , p . In this case, N denotes the dimension of the FEM space V h , p and equals the number of (interior with respect to Ω ) degrees of freedom. The operator A h in the finite element case is a map of V h , p into V h , p so that A h v h : = g h , where g h V h , p is the unique solution to
( g h , w h ) = A ( v h , w h ) , for   all w h V h , p .
Let us denote by { ϕ j } j = 1 N the standard nodal basis of V h , p . With respect to this basis, the operator A h is represented by the matrix
A = M 1 S , where S i , j = A ( ϕ i , ϕ j ) , M i , j = ( ϕ i , ϕ j ) .
In accordance with the terminology adopted in the finite element method, M and S are the mass and stiffness matrices, respectively. They are both symmetric and positive definite (SPD).
Obviously, if w = A h z and w , z R N are the coefficient vectors corresponding to w , z V h , p , then w = A z . Now, for the coefficient vector f h corresponding to f h = π h f , we obtain f h = M 1 F h , where F h is the vector with entries
F h , j = ( f , ϕ j ) , for j = 1 , 2 , , N .
Thus, using vector notation so that u h is the coefficient vector representing the solution u h through the nodal basis, we can write the finite element approximation of (3) in the form of a system
A α u h = f h = M 1 F h , which   implies S u h = F h .
Therefore, for the finite element approximation of the subdiffusion problem (7) we obtain the equation
M A α u h = F h , or   equivalently u h = A α M 1 F h .
We note that the matrix A R N × N is SPD with respect to the Euclidean space defined by the energy inner product associated with the mass matrix, that is
( v , w ) M = ( M v , w ) .
This follows directly from the definition (10), which implies the equality
( A v , w ) M = ( S v , w ) .
In the context of our study, the FEM problem (8) is by default large-scale, i.e., N > > 1 . When analyzing the new BURA-SD algorithm, we will assume that systems of the type ( S + d M ) v h = g h , d 0 and v h , g h R N can be solved numerically in an efficient way by a method that requires O ( N ) 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 O ( N ) as possible.
In the case of quasi-uniform partition T h and h FEM approximation, the condition number of the stiffness matrix S increases as h 2 when h 0 , i.e., κ ( S ) = O ( h 2 ) . Here, the focus is on the h p 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 S G L the standard nodal basis stiffness matrix and by S S P the stiffens matrix obtained if the corresponding spectral basis functions are used in the interior of the finite elements K T h p . If the mesh T h p is quasi-uniform, then the following estimates hold (see [29]):
κ ( S G L ) = O ( p 4 log p h 2 ) and κ ( S S P ) = O ( p 3 h 2 ) .
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 κ ( S G L ) when the number of layers increases. The behavior of S S P is significantly better. Under rather more general conditions, it was shown in [29] that, for meshes that are refined geometrically toward singularities, κ ( S S P ) is independent of the number of refinement levels.
The condition number of the mass matrix is also of interest when the BURACD solver is implemented. In the case of quasi-uniform partition T h and h FEM approximation, κ ( M ) = O ( 1 ) . Now, unlike (13), if T h p is quasi-uniform and the h p FEM is applied,
κ ( M G L ) = O ( p 4 log p ) and κ ( M S P ) = O ( p 3 ) .
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, κ ( M S P ) 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 κ ( A ) = κ ( M 1 S ) 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 κ ( A ) is needed as an input parameter.

3. BURA-SD

Let us denote by r α , k ( z ) the best uniform rational approximation (BURA) of degree k of z α in [ 0 , 1 ] , belonging to the diagonal of the Walsh table. This means that we consider the rational functions r k ( z ) = P k ( z ) / Q k ( z ) , where P k and Q k are polynomials of equal degree k. By definition, r α , k ( z ) is the error minimizer for which
E α , k : = min r k ( z ) R ( k , k ) max z [ 0 , 1 ] | z α r k ( z ) | .
The following additive and multiplicative representations of r α , k ( z ) can be used in the implementation of the BURA methods and algorithms:
r α , k ( z ) = c 0 + i = 1 k c i z d i = c 0 i = 1 k z ξ i z d i , where   c i > 0 .
As was shown in [13], all k zeros ξ 1 , , ξ k and poles d 1 , , d k are real and negative. Even more precisely, they are interlacing, that is, with appropriate numbering
0 > ξ 1 > d 1 > ξ 2 > d 2 > > ξ k > d k > .
The following asymptotically sharp exponential estimate of E α , k (see for more details [13] and references there in) holds:
E α , k 4 α + 1 sin ( α π ) e 2 π α k .
In accordance with the introduced notations, the BURA approximation of A α for α ( 0 , 1 ) is defined as
A α λ 1 , h α r α , k ( λ 1 , h A 1 ) ,
where λ 1 , h is the smallest eigenvalue of A . In the particular case of FEM approximation, see (10), λ 1 , h is the smallest positive eigenvalue of the generalized eigenvalue problem
S v = λ M v .
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 w h the BURA approximation of the solution u h of h p FEM system (11) defined as follows
u h : = λ 1 , h α r α , k ( λ 1 , h A 1 ) f h .
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 r α , k ( 1 / z ) , while the multiplicative one deals with the factorization of r α . k ( 1 / z ) . Thus, transforming the expressions (16), we obtain the equalities
r ˜ α , k ( z ) : = r α , k ( 1 / z ) = c ˜ 0 + i = 1 k c ˜ i z d ˜ i = c ˜ 0 i = 1 k z ξ ˜ i z d ˜ i .
Here, c ˜ i > 0 and d ˜ i , ξ ˜ i < 0 . In this way, the variants of BURA method based on (19) read as
u h = λ 1 , h α c ˜ 0 I + i = 1 k λ 1 , h c ˜ i ( A λ 1 , h d ˜ i I ) 1 f h , u h = λ 1 , h α c ˜ 0 i = 1 k ( A λ 1 , h ξ ˜ i I ) ( A λ 1 , h d ˜ i I ) 1 f h .
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 A . We should note here that the matrix A = M 1 S 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
u h = λ 1 , h α c ˜ 0 M 1 + i = 1 k λ 1 , h c ˜ i ( S λ 1 , h d ˜ i M ) 1 F h , u h = λ 1 , h α c ˜ 0 M 1 i = 1 k ( S λ 1 , h ξ ˜ i M ) ( S λ 1 , h d ˜ i M ) 1 F h .
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 S with the mass matrix M .
The discussed BURA method can be directly applied to the h p FEM discretization of (3) as follows:
  • BURA method
(1)
Consider the fractional diffusion problem (3), A α u = f , α ( 0 , 1 ) .
(2)
Choose the h p FEM space V h , p V and compute the stiffness matrix S corresponding to the elliptic operator A and the mass matrix M .
(3)
Choose the degree of rational approximation k and find the best uniform rational approximation r α , k ( z ) of z α , z [ 0 , 1 ] .
(4)
Applying (18), compute the BURA approximation of the solution u h of the system (11), A α u h = f h .
(5)
Find the BURA approximation of u V as the h p FEM function u h V h represented by the nodal basis vector u h .
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.
  • BURA-SD method
(1)
Consider the fractional diffusion problem (3), A α u = f , α ( 0 , 1 ) .
(2)
Choose the degree of rational approximation k and find the best uniform rational approximation r α , k ( z ) of z α , z [ 0 , 1 ] .
(3)
Find the semi-discrete BURA approximation u ^ of u by the formula
u ^ : = λ 1 α r α , k ( λ 1 A 1 ) f .
(4)
Choose the h p FEM space V h , p V and compute the stiffness matrix S corresponding to the elliptic operator A and the mass matrix M . Apply the h p FEM to (22) and obtain the approximation u ^ h of u ^ , which satisfies the equality
u ^ h : = λ 1 α r α , k ( λ 1 A 1 ) f h ,
(5)
Find the BURA approximation of u ^ V as the h p FEM function u ^ h V h represented by the nodal basis vector u ^ h .
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
lim λ 1 , h λ 1 u h = u ^ h , a n d   c o n s e q u e n t l y lim λ 1 , h λ 1 u h = u ^ h .
Let us recall that the first eigenvalue of the operator A associated with the elliptic boundary value problem (1) is uniformly separated from 0. Thus, λ 1 = O ( 1 ) and accordingly, for a sufficiently fine h p FEM discretization, λ 1 , h = O ( 1 ) . 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 E α , k 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.

4. Error Analysis

In this section, we estimate the accuracy of the numerical solution u ^ h obtained by BURA-SD method. The analysis follows the representation E = u u ^ h = ( u u ^ ) + ( u ^ u ^ h ) , which corresponds to the successive steps of the method.

4.1. Error of the Semi-Discrete BURA-SD Approximation

We begin by estimating the error of the semi-discrete BURA-SD approximation E 1 = u u ^ . The orthonormal basis of eigenfunctions of A is used to expand the right-hand side f, thus obtaining the following estimate of E 1 in the energy norm induced by A
| | E 1 | | A 2 = j = 1 λ j α ( f , ϕ j ) ϕ j λ 1 α r α , k ( λ 1 A 1 ) j = 1 λ j ( f , ϕ j ) ϕ j , j = 1 λ j α ( f , ϕ j ) ϕ j λ 1 α r α , k ( λ 1 A 1 ) j = 1 λ j ( f , ϕ j ) ϕ j A = λ 1 2 α j = 1 λ j λ j λ 1 α r α , k λ 1 λ j 2 ( f , ϕ j ) .
Then
| | E 1 | | A 2 = λ 1 , h 2 α j = 1 λ j θ j α r α , k ( θ j ) 2 ( f , ϕ j )
where θ j : = λ 1 / λ j ( 0 , 1 ] [ 0 , 1 ] . Now, we apply (15) to obtain
| | E 1 | | A 2 λ 1 2 α E α , k 2 j = 1 λ j ( f , ϕ j ) = λ 1 2 α E α , k 2 | | f | | A 2 .
The following lemma is an immediate consequence of this estimate.
Lemma 1.
For each α ( 0 , 1 ) the error of the semi-discrete BURA-SD approximation E 1 = u u ^ decreases exponentially as the degree of the best uniform rational approximation k increases, satisfying the estimate
| | E 1 | | H 1 c 1 λ 1 α 4 α + 1 sin ( α π ) e 2 π α k | | f | | H 1 .
Here, the constant c 1 > 0 is independent of α and k.
Proof. 
The upper bound in the inequality (26) follows from the equivalence of the norms | | . | | H 1 and | | . | | A , and from (25), where (17) is applied to estimate E α , k . □

4.2. h p FEM Approximation of Singularly Perturbed Diffusion-Reaction Equations

Without limiting the general applicability of the results, from now on, for convenience, we will analyze the additive representation of the BURA-SD method. So, let us consider the semi-discrete approximation in the form
u ^ = λ 1 α c ˜ 0 I + i = 1 k λ 1 c ˜ i ( A λ 1 d ˜ i I ) 1 f ,
where I is the identity operator, c ˜ i > 0 and d ˜ i < 0 are real constants. The above formula follows from (19) as a semi-discrete analogue of the first equality from (20). It is important for the analysis that the coefficients d ˜ i have an exponential growth when i tends to 1. This is explained by the exponential clustering toward 0 of the poles of r α , k .
In this context, we will refer to [26], where the behavior of d ˜ i is discussed, with extreme growth illustrated by numerical results. For example, there we see that for k = 45 the first nine coefficients d ˜ i for α = 1 / 4 are greater than 10 20 . For k = 70 the results are similar: for α = 1 / 4 and α = 1 / 2 , for the first 21 and the first 8, respectively, d ˜ i > 10 20 .
Now, we rewrite (27) in the form
u ^ = λ 1 α c ˜ 0 I + i = 1 k ( c ˜ i / d ˜ i ) ( ε i 2 A + I ) 1 f ,
where ε i 2 = 1 / ( λ 1 d ˜ i ) . Therefore, we will have to solve diffusion-reaction problems with small positive parameters ε i .
Thus, to ensure exponential convergence of the semi-discrete BURA-SD approximation, a crucial role is played by obtaining a robust error estimate when applying the h p FEM to singularly perturbed diffusion-reaction problems of the form
ε 2 · ( a ( x ) u ε ( x ) ) + u ε ( x ) = f ( x ) , for x Ω , u ε ( x ) = 0 , for x Ω ,
where ε > 0 is a small parameter. It is known (see, e.g., [35]) that u ε exhibits boundary layers near Ω with additional singularities induced by the corner points. A similar behavior is also typical for the solution of the fractional diffusion problem (3).
In the present study, we will use the same finite element space V h p for the numerical solution of all diffusion-reaction subproblems that appear in the sum in (27) and in (28), respectively. In constructing V h p , it is natural to apply both the geometric refinement of the boundary layer mesh and the geometric refinement of the corner mesh.
Here, we will use the meshes T h p = T g e o , σ L , n defined in [20]. They are generated by geometric refinement of a small number of so-called reference patches associated with the regular initial macro-triangulation T H of (possibly curvilinear) quadrilaterals. For non-trivial patches with integers L n , we denote the number of layers of refinements towards an edge and towards a vertex, respectively. Omitting some details of the strict definitions and simplifying the notations, for a given parameter σ ( 0 , 1 ) the construction of T g e o , σ L , n is based on the following assumptions:
(A1)
Trivial patches are not further refined.
(A2)
Patches with L layers of geometric refinement towards an edge are characterized by a reference macroelement [ 0 , 1 ] 2 containing L rectangles that are determined by the nodes ( 0 , 0 ) , ( 0 , 1 ) , and ( 0 , σ i ) , ( 1 , σ i ) , i = 0 , 1 , L .
(A3)
Patches with n layers of geometric refinement towards a corner are characterized by a reference macroelement [ 0 , 1 ] 2 containing triangles that are determined by the nodes ( 0 , 0 ) , and ( 0 , σ i ) , ( σ i , 0 ) , ( σ i , σ i ) , i = 0 , 1 , , n .
(A4)
Patches with L layers of geometric refinement towards an edge and n layers of geometric refinement towards a corner that is vertex of the same edge are characterized by a reference macroelement [ 0 , 1 ] 2 with a tensor product of refinements of the types described in assumptions (A2) and (A3), respectively.
(A5)
Patches with L layers of geometric refinement towards two adjacent edges and n layers of geometric refinement towards a corner that is a common vertex of these edges are characterized by a reference macroelement [ 0 , 1 ] 2 with a tensor product of two refinements described in assumption (A2) and a refinement described in assumption (A3), respectively.
Three examples of patches characterized by the corresponding reference macroelements in [ 0 , 1 ] 2 are shown in Figure 1. The assumptions (A1–A5) actually define a catalogue of possible refinements. Additional figures illustrating such patches can be found in [20].
The following lemma will be used in the upcoming analysis.
Lemma 2.
Let Ω R 2 be a curvilinear polygon with J vertices and let the mesh T g e o , σ L , n satisfies the assumptions (A1–A5). Then, there exist constants C, γ > 0 , β [ 0 , 1 ) (depending only on the data a ( x ) , f, and on the macro-triangulation T H ), so if ε ( 0 , 1 ] and L satisfies the scale resolution condition
σ L ε ,
then
| | u ε u h ε | | H 1 C ε 2 p 9 σ ( 1 β ) n + e γ p .
Here, u h ε V h , p H 1 ( Ω ) is the h p FEM numerical solution of (27), where the finite element space V h , p is defined on T g e o , σ L , n . In addition, the dimension of V h , p is bounded as follows
N : = d i m ( V h , p ) C N J p 2 L 2 + n
where the constant C N > 0 can depend only on a ( x ) , f and T H .
Proof. 
The error estimate (31) and the inequality (32) follow from Proposition 2 of [20]. To this end, the analysis must be reproduced in terms of the definitions and notations introduced above. This result can also be derived directly from the earlier paper [35]. □

4.3. Exponential Convergence of the Fully-Discrete Approximation

Let us write the error of the fully discrete approximation as a sum of the following two component errors
E = u u ^ h = ( u u ^ ) + ( u ^ u ^ h )
or
E = E 1 + E 2 , where E 1 = u u ^ and E 2 = u ^ u ^ h .
The E 1 component was estimated in Lemma 1. For the second term, we apply (28) to obtain the representation
E 2 = u ^ u ^ h = i = 1 k ( c ˜ i / d ˜ i ) u ^ ε i u ^ h ε i ,
where
u ^ ε i = ε i 2 A + b ˜ i I 1 f
and
ε i 2 = 1 λ 1 d ˜ i .
In (33), u ^ h ε i V h , p is the h p FEM approximation of u ^ ε i obtained using the mesh T g e o , σ L , n . Now, under the assumptions of Lemma 2 we obtain the estimate
| | E 2 | | H 1 C i = 1 k ( c ˜ i / d ˜ i ) ε i 2 p 9 σ ( 1 β ) n + e γ p = C i = 1 k c ˜ i λ 1 p 9 σ ( 1 β ) n + e γ p .
Lemma 3.
Let Ω R 2 be a curvilinear polygon with J vertices and let the mesh T g e o , σ L , n be constructed in accordance with assumptions (A1–A5). Also, let L be large enough so that the scale resolution condition
σ L 1 / ( λ 1 d ˜ 1 )
is satisfied. Then, there exist constants C 2 , γ > 0 , β [ 0 , 1 ) (depending only on the data a ( x ) , f ( x ) , and on the macro-triangulation T H ) such that
| | E 2 | | H 1 C 2 c ˜ 1 k p 9 σ ( 1 β ) n + e γ p .
Proof. 
Analogous to { ξ i } i = 1 k , { d i } i = i k , the coefficients { c ˜ i } i = 1 k , { d ˜ i } i = 1 k satisfy the inequalities
> c ˜ 1 > c ˜ 2 > > c ˜ k > 0
and
< d ˜ 1 < d ˜ 2 < < d ˜ k < 0 .
Now, we use the equality
ε i 2 = d i / λ i = 1 / ( λ i d ˜ i ) ,
and taking into account that { d ˜ i } i = 1 k increases monotonically, we obtain from (35) that
( σ L ) 2 < 1 / ( λ 1 d ˜ 1 ) < 1 / ( λ 1 d ˜ i ) .
Therefore, the scale resolution condition (30) is satisfied for each ε i , i = 1 , , k . Finally, the estimate (35) follows from (34), noting that λ 1 = O ( 1 ) and that { c ˜ i } i = 1 k increases monotonically. □
The following theorem completes our error analysis.
Theorem 1.
Let Ω R 2 be a curvilinear polygon with mesh T g e o , σ L , n constructed according to assumptions (A1–A5). Also, let L be large enough so that the scale resolution condition
σ L 1 / ( λ 1 d ˜ 1 )
is satisfied. Then, there exist constants C 2 , γ > 0 , β [ 0 , 1 ) (depending only on the data a ( x ) , f ( x ) , and on the macro-triangulation T H ) such that
| | u u ^ h | | H 1 C 1 e 2 π α k + C 2 c ˜ 1 k p 9 σ ( 1 β ) n + e γ p .
Proof. 
The estimate (38) is a direct consequence of Lemmas 2 and 3. □

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 h p FEM exponential accuracy with respect to the mesh parameters ( n , p ) , 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
1 d ˜ 1 c d e c d e k
with positive constants c d and c d e that can only depend on the fractional power α .
The exponential decrease of 1 / d ˜ i with increasing i for k = 30 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 1 / d ˜ 1 , taking into account the global monotonicity of d ˜ i . The BRASIL software [33] is used to calculate d ˜ i , see also [38].
Lemma 4.
Let
L 1 ln ( 1 / σ ) ln 1 λ 1 c d + c d e k .
Then, the scale resolution condition (37) is satisfied.
Proof. 
We write (40) in the form
ln 1 λ 1 c d + c d e k L ln ( 1 / σ ) .
Then
1 λ 1 c d e c d e k ( 1 / σ ) L
and therefore
σ L λ 1 c d e c d e k .
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 T g e o , σ L , n be large enough to provide the necessary h p FEM accuracy to the ε i 2 -perturbed reaction-diffusion problems that arise in the BURA-SD method. It follows from Lemma 4 that there exists a constant c k > that can only depend on α, so the scale resolution condition is satisfied if
L = c k k .
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
e 2 π α k k p 9 σ ( 1 + β ) n .
Therefore,
k ln k p + n ,
and consequently
k ( ln k p + n ) 2 .
Analogously, the sufficient condition for balancing the BURA-SD and p FEM errors has the form
e 2 π α k k p 9 e γ p .
Then,
k ln k p + p ,
and consequently
k ( ln k p + p ) 2 .
Using (43) and (44) we obtain the following lemma.
Lemma 5.
Let the scale resolution condition (37) be satisfied and let
k = n 2 + p 2 .
Then,
| | u u ^ h | | H 1 C e ( n 2 + p 2 ) p 9 σ ( 1 β ) n + e γ p .
Similar to Lemma 3, the constants C e , γ > 0 , β [ 0 , 1 ) depend only on the data a ( x ) ,   f ( x ) and on the macrotriangulation T H .
We conclude our analysis with the following theorem.
Theorem 2.
Let Ω R 2 be a curvilinear polygon with mesh T g e o , σ L , n constructed according to assumptions (A1–A5). Then, there exist a mesh refinement parameter σ ( 0 , 1 ) , BURA-SD degree k, and constants C e , C c , γ > 0 , β [ 0 , 1 ) (depending only on the data a ( x ) , f ( x ) , and on the macro-triangulation T H ), such that the following error estimate holds
| | u u ^ h | | H 1 C e ( n 2 + p 2 ) p 9 σ ( 1 β ) n + e γ p
and the computational complexity is bounded by
N BURA C c p 4 L ( L 2 + n ) .
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 S λ 1 , h d ˜ i M R N × N , i = 1 , 2 , , k . 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 n i t MG = O ( 1 ) is uniformly bounded regardless of p. Then, the following estimate of the computational complexity N BURA of the BURA-SD method holds
N BURA = O k p 2 N ,
where p 2 represents the number of nonzero off-diagonal elements of the matrices. Thus, from (32) we obtain
N BURA = O k p 4 ( L 2 + n ) ,
and finally applying (42) we obtain
N BURA = O p 4 L ( L 2 + n ) .

6. Concluding Remarks

In this paper, we have developed a new method for the numerical solution of sub-diffusion problems in curvilinear polygons Ω R 2 . This means that there are practically no restrictions on the geometry of the bounded computational domain. The spectral definition of the elliptic operator A α , α ( 0 , 1 ] , is used. The proposed methodology allows the construction of a method that fully utilizes the strengths of the h p 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 O ( σ ( 1 β ) n + e γ p ) of the h p 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 N BURA = O ( p 2 L N ) . We recall that in constructing the h p 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 Ω × ( 0 , ) 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 M decoupled singularly perturbed diffusion-reaction equations in Ω . The h p 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 A α . In the context of the present paper, the integer M corresponds to the BURA-SD degree k.
In this work, we elaborated a novel approach to the BURA (here called BURA-SD) approximation of A α that is adapted to the specifics in the analysis of the h p 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 M 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 A α 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 h p 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 κ ( A ) .
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 O ( k ) . 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 h p FEM discretization case. In perspective, however, we are analyzing alternative possibilities to further develop our research in the spirit of meshless methods.

Funding

We acknowledge the support of the Centre of Excellence in Informatics and ICT, Grant No BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program and co-financed by the European Union through the European Structural and Investment funds. The work is also partially supported by the Bulgarian National Science Fund under grant No. KP-06-N72/2.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The author declares no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
BURABest Uniform Rational Approximation
BURA-SDBURA-Semi-Discrete
BURA-ARBURA-Additive Reduced
BURA-MRBURA-Multiplicative Reduced
FEMFinite Element Method
AP#Approach#
A#Assumption#

References

  1. Caffarelli, L.; Silvestre, L. An extension problem related to the fractional Laplacian. Comm. Partial. Differ. Equ. 2007, 32, 1245–1260. [Google Scholar] [CrossRef]
  2. Caffarelli, L.A.; Stinga, P.R. Fractional elliptic equations, Caccioppoli estimates and regularity. Ann. Inst. H. Poincaré C Anal. Non Linéaire 2016, 33, 767–807. [Google Scholar] [CrossRef]
  3. Chen, L.; Nochetto, R.H.; Otárola, E.; Salgado, A.J. Multilevel methods for nonuniformly elliptic operators and fractional diffusion. Math. Comp. 2016, 85, 2583–2607. [Google Scholar] [CrossRef]
  4. Kosturski, N.; Margenov, S. Analysis of BURA and BURA-based approximations of fractional powers of sparse SPD matrices. Fract. Calc. Appl. Anal. 2024, 27, 706–724. [Google Scholar] [CrossRef]
  5. Nochetto, R.H.; Otárola, E.; Salgado, A.J. A PDE approach to space-time fractional parabolic problems. SIAM J. Numer. Anal. 2016, 54, 848–873. [Google Scholar] [CrossRef]
  6. Bonito, A.; Borthagaray, J.P.; Nochetto, R.H.; Otárola, E.; Salgado, A.J. Numerical methods for fractional diffusion. Comput. Vis. Sci. 2018, 19, 19–46. [Google Scholar] [CrossRef]
  7. Vabishchevich, P.N. Numerically solving an equation for fractional powers of elliptic operators. J. Comput. Phys. 2015, 282, 289–302. [Google Scholar] [CrossRef]
  8. Vabishchevich, P.N. Numerical Solution of Nonstationary Problems for a Space-Fractional Diffusion Equation. Fract. Calc. Appl. Anal. 2016, 19, 116–139. [Google Scholar] [CrossRef]
  9. Vabishchevich, P.N. Approximation of a fractional power of an elliptic operator. Numer. Linear Algebra Appl. 2020, 27, e2287. [Google Scholar] [CrossRef]
  10. Bonito, A.; Lei, W.; Pasciak, J.E. Numerical approximation of the integral fractional Laplacian. Numer. Math. 2019, 142, 235–278. [Google Scholar] [CrossRef]
  11. Bonito, A.; Pasciak, J.E. Numerical approximation of fractional powers of elliptic operators. Math. Comp. 2015, 84, 2083–2110. [Google Scholar] [CrossRef]
  12. Harizanov, S.; Lazarov, R.; Margenov, S.; Marinov, P.; Pasciak, J. Analysis of numerical methods for spectral fractional elliptic equations based on the best uniform rational approximation. J. Comput. Phys. 2020, 408, 109285. [Google Scholar] [CrossRef]
  13. Harizanov, S.; Lazarov, R.; Margenov, S.; Marinov, P.; Vutov, Y. Optimal solvers for linear systems with fractional powers of sparse SPD matrices. Numer. Linear Algebra Appl. 2018, 25, e2167. [Google Scholar] [CrossRef]
  14. Hofreither, C. A unified view of some numerical methods for fractional diffusion. Comput. Math. Appl. 2020, 80, 332–350. [Google Scholar] [CrossRef]
  15. Harizanov, S.; Lazarov, R.; Margenov, S. A survey on numerical methods for spectral space-fractional diffusion problems. Fract. Calc. Appl. Anal. 2020, 23, 1605–1646. [Google Scholar] [CrossRef]
  16. Aceto, L.; Novati, P. Fast and accurate approximations to fractional powers of operators. IMA J. Numer. Anal. 2022, 42, 1598–1622. [Google Scholar] [CrossRef]
  17. Aceto, L.; Novati, P. Exponentially convergent trapezoidal rules to approximate fractional powers of operators. J. Sci. Comput. 2022, 91, 55. [Google Scholar] [CrossRef]
  18. Barakitis, N.; Ekström, S.E.; Vassalos, P. Preconditioners for fractional diffusion equations based on the spectral symbol. Numer. Linear Algebra Appl. 2022, 29, e2441. [Google Scholar] [CrossRef]
  19. Bertaccini, D.; Durastante, F. Nonlocal diffusion of variable order on complex networks. Int. J. Comput. Math. Comput. Syst. Theory 2022, 7, 172–191. [Google Scholar] [CrossRef]
  20. Banjai, L.; Melenk, J.M.; Schwab, C. Exponential convergence of hp FEM for spectral fractional diffusion in polygons. Numer. Math. 2023, 153, 1–47. [Google Scholar] [CrossRef]
  21. Duan, B. Padé-parametric FEM approximation for fractional powers of elliptic operators on manifolds. arXiv 2022. [Google Scholar] [CrossRef]
  22. Danczul, T.; Hofreither, C. On rational Krylov and reduced basis methods for fractional diffusion. J. Numer. Math. 2022, 30, 121–140. [Google Scholar] [CrossRef]
  23. Danczul, T.; Schöberl, J. A reduced basis method for fractional diffusion operators I. Numer. Math. 2022, 151, 369–404. [Google Scholar] [CrossRef]
  24. Denich, E.; Novati, P. A Gaussian Method for the Square Root of Accretive Operators. Comput. Methods Appl. Math. 2023, 23, 127–143. [Google Scholar] [CrossRef]
  25. Duan, B.; Zongze, Y. A quadrature scheme for steady-state diffusion equations involving fractional power of regularly accretive operator. SIAM J. Sci. Comput. 2023, 45, A2226–A2249. [Google Scholar] [CrossRef]
  26. Harizanov, S.; Kosturski, N.; Lirkov, I.; Margenov, S.; Vutov, Y. Reduced Multiplicative (BURA-MR) and Additive (BURA-AR) Best Uniform Rational Approximation Methods and Algorithms for Fractional Elliptic Equations. Fractal Fract. 2021, 5, 61. [Google Scholar] [CrossRef]
  27. Harizanov, S.; Lirkov, I.; Margenov, S. Rational Approximations in Robust Preconditioning of Multiphysics Problems. Mathematics 2022, 10, 780. [Google Scholar] [CrossRef]
  28. Harizanov, S.; Lazarov, R.; Margenov, S.; Marinov, P. Numerical solution of fractional diffusion-reaction problems based on BURA. Comput. Math. Appl. 2020, 80, 316–331. [Google Scholar] [CrossRef]
  29. Melenk, J.M. On condition numbers in hp-FEM with Gauss-Lobatto-based shape functions. J. Comput. Appl. Math. 2002, 139, 21–48. [Google Scholar] [CrossRef]
  30. Gunatilake, J. An improved multigrid solver for the p-hierarchical basis finite element method using a space decomposition smoother. Comput. Math. Appl. 2022, 124, 52–62. [Google Scholar] [CrossRef]
  31. Vejchodský, T.; Šolín, P. Improving conditioning of hp-FEM. In SNA’07 Modelling and Simulation of Challenging Engineering Problems; Institute of Geonics Czech Academy of Sciences: Ostrava, Czech Republic, 2007; pp. 126–129. [Google Scholar]
  32. Vejchodský, T.; Šolín, P. Static condensation, partial orthogonalization of basis functions, and ILU preconditioning in the hp-FEM. J. Comput. Appl. Math. 2008, 218, 192–200. [Google Scholar] [CrossRef]
  33. Hofreither, C. Software BRASIL. Available online: https://www.softwaresbrasil.com(accessed on 1 July 2024).
  34. Li, Y.; Zikatanov, L.; Zuo, C. A Reduced Conjugate Gradient Basis Method for Fractional Diffusion. SIAM J. Sci. Comput. 2024, S68–S87. [Google Scholar] [CrossRef]
  35. Melenk, J.M.; Schwab, C. hp FEM for reaction-diffusion equations. I. Robust exponential convergence. SIAM J. Numer. Anal. 1998, 35, 1520–1557. [Google Scholar] [CrossRef]
  36. Stahl, H.R. Best uniform rational approximation of xα on [0,1]. Acta Math. 2003, 190, 241–306. [Google Scholar] [CrossRef]
  37. Trefethen, L.N.; Nakatsukasa, Y.; Weideman, J.A.C. Exponential node clustering at singularities for rational approximation, quadrature, and PDEs. Numer. Math. 2021, 147, 227–254. [Google Scholar] [CrossRef] [PubMed]
  38. Hofreither, C. An algorithm for best rational approximation based on barycentric rational interpolation. Numer. Algorithms 2021, 88, 365–388. [Google Scholar] [CrossRef]
  39. Bulle, R.; Barrera, O.; Bordas, S.P.; Chouly, F.; Hale, J.S. An a posteriori error estimator for the spectral fractional power of the Laplacian. Comp. Meth. Appl. Mech. Eng. 2023, 407, 115943. [Google Scholar] [CrossRef]
Figure 1. Patches with: (left) L = 4 layers of geometric refinement towards an edge (A2); (centre) n = 4 layers of geometric refinement towards a corner (A3); (right) L = 4 layers of geometric refinement towards an edge and n = 4 layers of geometric refinement towards a corner that is vertex of the same edge (A4).
Figure 1. Patches with: (left) L = 4 layers of geometric refinement towards an edge (A2); (centre) n = 4 layers of geometric refinement towards a corner (A3); (right) L = 4 layers of geometric refinement towards an edge and n = 4 layers of geometric refinement towards a corner that is vertex of the same edge (A4).
Mathematics 12 02266 g001
Figure 2. Behavior of ln 1 d ˜ i for 1 i { 1 , 2 , , 30 } and α { 0.25 , 0.50 , 0.75 } .
Figure 2. Behavior of ln 1 d ˜ i for 1 i { 1 , 2 , , 30 } and α { 0.25 , 0.50 , 0.75 } .
Mathematics 12 02266 g002
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Margenov, S. Exponential Convergence and Computational Efficiency of BURA-SD Method for Fractional Diffusion Equations in Polygons. Mathematics 2024, 12, 2266. https://doi.org/10.3390/math12142266

AMA Style

Margenov S. Exponential Convergence and Computational Efficiency of BURA-SD Method for Fractional Diffusion Equations in Polygons. Mathematics. 2024; 12(14):2266. https://doi.org/10.3390/math12142266

Chicago/Turabian Style

Margenov, Svetozar. 2024. "Exponential Convergence and Computational Efficiency of BURA-SD Method for Fractional Diffusion Equations in Polygons" Mathematics 12, no. 14: 2266. https://doi.org/10.3390/math12142266

APA Style

Margenov, S. (2024). Exponential Convergence and Computational Efficiency of BURA-SD Method for Fractional Diffusion Equations in Polygons. Mathematics, 12(14), 2266. https://doi.org/10.3390/math12142266

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop