Next Article in Journal
A Formal Approach to Coercion Resistance and Its Application to E-Voting
Next Article in Special Issue
On a Framework for the Stability and Convergence Analysis of Discrete Schemes for Nonstationary Nonlocal Problems of Parabolic Type
Previous Article in Journal
Modified Mann Subgradient-like Extragradient Rules for Variational Inequalities and Common Fixed Points Involving Asymptotically Nonexpansive Mappings
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Rational Approximations in Robust Preconditioning of Multiphysics Problems

Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(5), 780; https://doi.org/10.3390/math10050780
Submission received: 31 January 2022 / Revised: 22 February 2022 / Accepted: 25 February 2022 / Published: 28 February 2022

Abstract

:
Multiphysics or multiscale problems naturally involve coupling at interfaces which are manifolds of lower dimensions. The block-diagonal preconditioning of the related saddle-point systems is among the most efficient approaches for numerically solving large-scale problems in this class. At the operator level, the interface blocks of the preconditioners are fractional Laplacians. At the discrete level, we propose to replace the inverse of the fractional Laplacian with its best uniform rational approximation (BURA). The goal of the paper is to develop a unified framework for analysis of the new class of preconditioned iterative methods. As a final result, we prove that the proposed preconditioners have optimal computational complexity O ( N ) , where N is the number of unknowns (degrees of freedom) of the coupled discrete problem. The main theoretical contribution is the condition number estimates of the BURA-based preconditioners. It is important to note that the obtained estimates are completely analogous for both positive and negative fractional powers. At the end, the analysis of the behavior of the relative condition numbers is aimed at characterizing the practical requirements for minimal BURA orders for the considered Darcy–Stokes and 3D–1D examples of coupled problems.

1. Introduction

This study is devoted to the development of computationally efficient methods for numerically solving multiphysics and multiscale problems. A representative set of such examples is discussed in the review article [1]. Such applications are also considered in [2,3,4]. More specifically, our attention is focused on a class of problems that involve coupling at interfaces, which are manifolds of lower dimensions. Such mathematical models are often described by saddle-point systems of differential equations. For example, this is the case when either a mixed formulation is used to ensure better accuracy and robustness or when Lagrangian multipliers are introduced to interpret the coupling at the interface. Thus, interface conditions in fractional Sobolev spaces naturally appear. For more related details, we can refer to [5,6,7].
Aimed at presumably more realistic applications, we are interested in the numerical solutions of multidimensional boundary value problems in bounded domain with general geometry Ω R d . We assume that the geometry of the interface Γ Ω is general, as well. After discretization by a certain mesh method (such as finite differences, finite elements, finite volumes, etc.) on a properly generated unstructured mesh (with possible mesh refinements), we obtain a system of linear algebraic equations with a sparse unstructured matrix A . Under the above assumptions, the discrete problems are apparently very large scale. For such problems, the advantages of iterative solution methods (solvers) are well expressed. Note that the condition number of the corresponding matrix can be very large. It could be due not only to the number of degrees of freedom but also to the structure of the mesh. For the multiscale and multiphysical problems under consideration, the media can be very highly heterogeneous and the anisotropy ratio can be very large. Thus, looking for solvers of optimal computational complexity gives rise to the requirement for efficient and robust preconditioning.
The block-diagonal preconditioning is among the most efficient approaches for numerically solving the studied saddle-point problems. Let us consider the system of differential equations A U = F written in the block-operator form
A 1 , 1 A 1 , 2 A 1 , n Ω A 1 , Γ A 2 , 1 A 2 , 2 A 2 , n Ω A 2 , Γ A n Ω , 1 A n Ω , 2 A n Ω , n Ω A n Ω , Γ A Γ , 1 A Γ , 2 A Γ , n Ω A Γ , Γ U = F ,
where the first n Ω rows correspond to the volume equations while the last one represents the interface condition.
Definition 1.
The block-diagonal preconditioners have zero off-diagonal blocks and nonsingular diagonal blocks. They are also known as additive preconditioners.
In this study, we analyze two types of block preconditioners that have optimal convergence rate, that is, the number of iterations is O ( 1 ) . Following the notation from Equation (1), they can be written as follows:
C ( 1 ) = d i a g A 1 ( 1 ) , A 2 ( 1 ) , , A n Ω ( 1 ) , ( Δ Γ + I Γ ) α ,
and
C ( 2 ) = d i a g A 1 ( 2 ) , A 2 ( 2 ) , , A n Ω ( 2 ) , I Γ , ( Δ Γ ) β ,
where α ( 0 , 1 ) , respectively β ( 1 , 0 ) . The first n Ω blocks of C ( 1 ) and C ( 2 ) are selfadjoint positive definite differential operators, I Γ is the identity operator at the interface Γ , and Δ Γ stands for the Laplacian at Γ .
At discrete level, the implementation of the first n Ω blocks of C ( 1 ) requires solving systems with sparse symmetric and positive definite (SPD) matrices. The same holds true for the first n Ω + 1 blocks of C ( 2 ) . For all of them, a standard PCG (e.g., AMG) solver of optimal computational complexity can be applied. The matrices, corresponding to the last blocks, are dense.
Thus, for the application of the preconditioner C ( 1 ) , the solution of a system with a matrix in the form ( A Δ Γ + I ) α is required, while the implementation of C ( 2 ) leads to matrix vector multiplication with A Δ Γ β . Here, A Δ Γ stands for the matrix corresponding to the discretization of Δ Γ . In both cases, numerical computations with a fractional power in the unit interval ( 0 , 1 ) of an SPD matrix are involved—either via solving a linear system of equations or via a matrix-vector multiplication. As we will see later on, although the implementation of the last block of C ( 2 ) is inverse free, it is not less challenging (but even more complicated) than the case of C ( 1 ) .
In this context, multigrid methods for discrete fractional Sobolev spaces are developed in [5]. The abstract additive multilevel framework is adapted there. The constructed smoother involves small ( ν + 1 ) × ( ν + 1 ) blocks in the form A k , ν α , where ν is the graph-degree of the mesh node k.
As an alternative, we propose to use the best uniform rational approximation (BURA) method for A α [8,9,10,11]. The aim is to develop a robust preconditioning algorithm of optimal computational complexity. As correctly noted in [5], a disadvantage of the method introduced in 2018 in [11] was that the accuracy depended on the condition number of the matrix A . This drawback has been overcome in the improved BURA approximation developed in [10]; see also the survey paper [8]. It reads as A α λ 1 , h α r α , k ( λ 1 , h A 1 ) , where r α , k is the best uniform rational approximation of degree k of z α on [ 0 , 1 ] and λ 1 , h is the smallest eigenvalue of A . The L 2 error of the BURA approximation decreases with k, satisfying the asymptotic estimate E α , k = O e 2 π k α . Further improvement of the computational efficiency of the method is presented in [9], where problem-specific model reduction techniques are utilized.
The structure of the paper is as follows. Two examples of coupled multiphysics problems, leading to the class of discussed block-diagonal preconditioners, are presented in Section 2. Section 3 is devoted to the construction and the basic properties of the BURA method. Then, the accuracy of the BURA-based preconditioners is analyzed in Section 4, including the derived estimates for the corresponding condition number. The theoretical results are supported by representative numerical experiments. The implementation issues are discussed in Section 5, ending with the proof of optimal computational complexity of the studied methods and algorithms. The motivation for the choice of numerical tests presented in Section 6 is twofold. On one hand, they illustrate the accuracy of the derived theoretical estimates, while on the other, they provide an additional practical insight for the efficient usage of BURA preconditioners with a small degree k. Finally, we give some conclusions in Section 7.

2. Monolithically Coupled Multiphysics Problems: Two Examples

In this section, we present two examples of related multiphysics problems that lead to block-diagonal preconditioners in the form (2) and (3). These preconditioners are already known and they are important for the analysis, presented in this paper as real world examples, where fractional elliptic problems are formulated. In this paper, we do not analyze the specific properties of the approximations of multiphysics equations or the influence of the different blocks of the preconditioner on the total condition number of the preconditioned system. The main goal of our study is to obtain detailed and asymptotically exact estimates of the accuracy of the BURA method. To simplify notation, throughout the paper we will use N for both the size of an SPD matrix and the number of unknowns of a coupled discrete problem. From the context, it will be clear to which of the two possibilities N refers to. Similarly, we denote by A various SPD matrices and again the context will clarify the particular object of interest.

2.1. Example 1: Coupling of Darcy and Stokes Flows

Let Ω be a closed domain in R 2 , which is a disjoint union of two subdomains Ω S and Ω D , that is, Ω S Ω D = . Figure 1 shows two examples of model problems with simple geometry and more realistic problem with a more general geometry of the interface.
The exchange through the interface Γ between a Stokes flow in Ω S and a Darcy flow in Ω D is modeled by the system of equations in the form [7,12,13]:
μ Δ u S + p S = f S in Ω S , · u S = 0 , in Ω S , K 1 u D + p D = 0 in Ω D , · u D = f D in Ω D , u D · n u S · n = 0 on Γ , μ u S n · n + p S = p D on Γ , μ u S n · t D u S · t = 0 on Γ .
Here, u S , p S are the velocities and pressure for the Stokes problem in Ω S , u D , p D are the velocities and pressure for the Darcy problem in Ω D , K and μ are the hydraulic conductivity, respectively, the fluid viscosity, while the coefficient D corresponds to the Beavers–Joseph–Saffman condition. Finally, n and t are unit vectors that are normal and tangential to the interface Γ . We further assume that the problem is equipped with correct boundary conditions.
Similarly to (1), the system of differential Equation (4) can be written in the block-operator form A D S U = F , where
A D S = A 1 , 1 A 1 , 2 0 0 0 0 0 A 2 , 1 0 0 0 0 0 0 0 0 A 3 , 3 A 3 , 4 0 0 0 0 0 A 4 , 3 0 0 0 0 A 5 , 1 0 A 5 , 3 0 A 5 , 5 0 0 A 6 , 1 A 6 , 2 0 A 6 , 4 A 6 , 5 A 6 , 6 A 6 , 7 A 7 , 1 0 0 0 A 7 , 5 0 0 .
The block sparsity structure is shown here. Then, after discretization by some mesh method (such as finite elements, finite differences, or finite volumes) we will obtain a linear system, where all blocks that correspond to nonzero operators in (5) will be also sparse.
There are two groups of methods for preconditioned iterative solution of coupled multiphysics problems: monolithic and based on some splitting [5,14,15]. For example, splitting techniques include operational splitting (by physical processes), domain decomposition (geometrical splitting), and block partitioning (algebraic splitting). The methods considered in this paper belong to the second group, following the motivation from [13]. Additionally, our study is supported by some recent achievements in the solution methods for fractional diffusion problems, or more generally for linear algebraic systems with fractional power of sparse SPD matrices.
The case μ = 1 is considered in [13]. It should be noted that this simplified problem is representative enough to illustrate why the norms studied in [12] are not sufficient for the proposed preconditioning approach. In [16], it is shown that the block-diagonal preconditioner
C D S = Δ + D T t T T t K 1 ( · + I ) I K I ( Δ + I ) 1 / 2 ,
provides an optimal condition number estimate. Here, T t is a tangential trace operator, while I stands for the identity operator. The first four blocks of (6) are standard components for preconditioning Darcy and Stokes problems. Thus, after discretization, the matrices corresponding to the first two blocks are sparse SPD, to which some Preconditioned Conjugate Gradients (PCG) solver (e.g., AMG) of optimal computational complexity can be applied. Note that the local operator T t T T t only increases the diagonal dominance in the first block. The last block leads to a system with the dense matrix ( A Δ Γ + I ) 1 / 2 .
This example raises the question of the efficient solution/preconditioning of systems with matrices in the form A α , A is sparse SPD, α ( 0 , 1 ) . Until recently, this was a challenge.

2.2. Example 2: Coupling a 1D Structure Embedded in Ω R 3

Let Ω R 3 be a bounded domain, while Γ represents a 1D manifold (structure) inside Ω , see Figure 2. The following trace coupled problem is considered [17]:
Δ u + u + p δ Γ = f in Ω , Δ v + v p = g on Γ , T u v = h on Γ .
Here u, v, p are unknowns, the 1D manifold Γ is parameterized in terms of s, the Laplacian on Γ is understood in terms of the Laplace–Beltrami operator, T : Ω Γ is a suitable trace operator, and p δ Γ is a Dirac measure such that Ω p ( x ) δ Γ w ( x ) d x = Γ p ( s ) δ Γ w ( s ) d s for a continuous function w.
At the same time, it is worth noting that the equations used in the real-life simulations are more complicated than (7); thus, the above example can be considered as a relevant model problem.
As it can be seen in [17], the standard finite element methods provide accurate error estimates in the norms of the introduced fractional order Sobolev spaces. The saddle-point matrix corresponding to the obtained system of linear algebraic equations is sparse. Arguments, similar to those discussed after (4) in the case of the Darcy–Stokes problem, are applicable here. However, standard techniques are not directly applicable to effective preconditioning. In what follows, we will consider the block-diagonal preconditioner
C 3 1 = Δ + I Δ + I Δ 0.14 ,
introduced in [6]. The choice of the fractional power in the third block of (8) is analyzed there, observing that values α [ 0.145 , 0.1 ] yield bounded condition numbers. The presented additional stability analysis has led to the preferred value α = 0.14 . Similar arguments to those concerning the implementation of the preconditioner C D S are applicable to C 3 1 .
The efficient approximate matrix vector multiplication with the matrix A α , where A is a sparse SPD matrix and α ( 0 , 1 ) , is the challenging question raised by this example.

2.3. Implementation of C D S and C 3 1 : The Key Question

The key questions in implementing the block-diagonal preconditioners C D S and C 3 1 are raised by the last blocks. They concern the fractional powers of SPD matrices defined through the spectral decomposition
A = W T D W ,
where the columns of W are the orthonormal eigenvectors of A and D stands for the diagonal matrix of the related eigenvalues. Under these assumptions, A W T = W T D , the matrix W is unitary, W T W = I , which leads directly to (9). Therefore
A α = W T D α W , A α = W T D α W .
The Formula (10) could be used in practical computations if the N eigenvectors and eigenvalues are explicitly known and if the matrix vector multiplication with W is equivalent to a fast Fourier transform. In such cases, the computational complexity is almost linear, O ( log N ) . However, when talking about fractional Laplacian on Γ , this limits the applications to interfaces with simple geometry (like the left sides of Figure 1 and Figure 2), uniform meshes, and to the lowest order finite element approximations.
The multigrid methods proposed and analyzed in [5] are aimed at handling general geometries of the interface Γ . As noted in the Introduction, the smoother includes non-local (dense) but small ( ν + 1 ) × ( ν + 1 ) blocks A k , ν α , where ν is the graph-degree of the mesh node k. Under certain regularity assumptions and fixed number of the multigrid levels J, the condition number κ is uniformly bounded for each given α ( 0 , 1 ) . More precisely, κ = O ( log N ) or κ = O ( log 2 N ) , depending on the regularity assumptions.
The case of negative power β ( 1 , 0 ) corresponds to the preconditioner (3), and more specifically to (8). The implementation of this kind of preconditioner requires a matrix vector multiplication with A β . For this multiplication, one can use the representation
A β = A A α , α = β + 1 ,
thus reducing the problem to numerically solving linear system with A α , α ( 0 , 1 ) . Following (11), the second variant of the multigrid method in [5] is adapted to handle negative powers. The condition number of the proposed preconditioner depends on N, including the case of a fixed number of multigrid levels J.
In the following sections, we consider the application of the BURA method in the general framework of the preconditioners C ( 1 ) and C ( 2 ) (Equations (2) and (3)) and, respectively, C D S and C 3 1 (Equations (6) and (8)). Our study is aimed at developing uniformly bounded condition number estimates, which are robust in terms of the geometry of the Γ interface, as well as in terms of regularity of the solution. The ultimate goal is to obtain preconditioned solution methods of optimal computational complexity.

3. The BURA Method

The fractional diffusion problems are non-local and their numerical solution is costly. The survey paper [8] tracks some of the last years’ achievements on this topic. They include a variety of methods based on: (i) extension of the problem from Ω R d to a local elliptic problem in Ω × ( 0 , ) , (ii) reformulation as a pseudo-parabolic problem in the cylinder Ω × ( 0 , 1 ) , (iii) methods based on approximation of the Dunford–Taylor integral representation of the solution, and (iv) methods based on the best uniform rational approximation (BURA) of z α on [ 0 , 1 ] . Applying rational approximations is one of the most successful approaches to efficiently solve non-local problems. For example, following the reformulation to a pseudo-parabolic problem, an efficient method based on diagonal Padé approximation to ( 1 + z ) α is developed in [18]. Recently, a unified view to the methods based on the approaches (i–iv) was presented in [19], showing that all of them can be interpreted as generated by some rational approximation of z α on [ 0 , 1 ] . In this sense, BURA methods are asymptotically the best.
The BURA method is originally introduced in [11] where the solution of the linear algebraic system
A α u = f , 0 < α < 1 ,
is written in the form u = A s α F with F = A s f , with s a (presumably small) positive integer. The proposed approximate solution method is based on the best uniform rational approximation (BURA) of the function z s α for 0 z 1 . Quite promising numerical results for s = 1 , accompanied by representative comparison analysis with other known methods, were presented in [11]. However, the convergence rate of the first BURA method was not robust with respect to the condition number of A , which is correctly noted in [5]. This shortcoming is completely overcome in [10]. Several aspects of the step-by-step development of the BURA method can be traced in the survey paper [8]. Note that in [19], it was shown that the available numerical solvers of (12) are based on uniform rational approximation of z α , z [ 0 , 1 ] . Therefore, methods based on best uniform rational approximation are expected to be the most efficient. From now on, when we talk about BURA, we will have in mind the following variant of the method proposed and analyzed in [10].
Let r α , k ( z ) be the best uniform rational approximation of z α in [ 0 , 1 ] in the class of rational functions r k ( z ) = P k ( z ) / Q k ( z ) , P k and Q k are polynomials of degree k. The case of rational approximations with equal degrees of the polynomials P and Q is known as the diagonal sequence of the Walsh table. The minimizer r α , k ( z ) satisfies the following sharp exponential error estimate with respect to the degree k [20]:
E α , k : = min r k ( z ) R ( k , k ) max z [ 0 , 1 ] | z α r k ( z ) | 4 α + 1 sin ( α π ) e 2 π α k .
Then, the BURA approximation of A α is defined as
A α λ 1 , h α r α , k ( λ 1 , h A 1 ) ,
where λ 1 , h is the smallest eigenvalue of the SPD matrix A .
The Padé approximation is also defined in the class of rational functions. However, while the minimizer r α , k ( z ) is implicitly defined, the Padé approximation is a Hermite interpolation. The advantage of the Padé approximation is that it is easier to compute.
Let A be obtained after finite difference discretization of a given elliptic operator. In the context of this paper, we can assume that A corresponds to the Laplacian with homogeneous Dirichlet boundary conditions. Then, λ 1 , h = O ( 1 ) , which follows, e.g., from the Poincaré–Steklov inequality. The same is true if the finite element method (FEM) is applied to discretize the Laplacian. In this case, A = M 1 K where K and M are the stiffens and mass matrices, respectively. Then, the SPD property of the FEM matrix A is understood in terms of the inner product induced by the mass matrix M .
Let us denote by
u k = λ 1 , h α r α , k ( λ 1 , h A 1 ) f
the BURA approximation of the solution u of the fractional diffusion linear system (12). Then, applying spectral basis representation of u , u k and of f , and the Parseval’s formula, we obtain the error estimate (see, for more details [10])
| | u u k | | 2 λ 1 , h α E α , k | | f | | 2 λ 1 , h α 4 α + 1 sin ( α π ) e 2 π α k | | f | | 2 .
The estimate (14) clearly shows the exponential rate with respect to the degree k of the rational approximation. Moreover, its dependence on the fractional power is explicit. Smaller α means lower accuracy. Since λ 1 , h = O ( 1 ) is uniformly bounded from below, the estimate is robust with respect to the condition number κ ( A ) . The estimate (14) is very sharp. The numerical tests conducted in [9] illustrate how close the computed errors are to the theoretical estimate, when the latter is greater than the double-precision accuracy.
The properties of the single roots and the related interlacing poles of r α , k ( z ) are well studied [21]. Both the additive and the multiplicative representation of the rational function r α , k ( 1 / z ) can be used in the implementation of the BURA method. The corresponding matrix valued expressions read as
u k = λ 1 , h α c ˜ 0 I + i = 1 k ( λ 1 , h c ˜ i ) ( A λ 1 , h d ˜ i I ) 1 f ,
and
u k = λ 1 , h α c ˜ 0 i = 1 k ( A λ 1 , h ξ ˜ i I ) ( A λ 1 , h d ˜ i I ) 1 f ,
where c ˜ i > 0 and d ˜ i , ξ ˜ i < 0 . A numerically stable approach for the computation of these coefficients is developed in [22]. Both the additive (15) and the multiplicative (16) algorithms reduce the problem to solving k systems with sparse SPD matrices A λ 1 , h d ˜ i I , which are positive diagonal perturbations of A . This leads to the understanding of the degree k as a universal measure for algorithmic efficiency, thus giving rise to the following asymptotic estimate of computational efficiency of the BURA methods
N α , k B U R A = O ( k N ) .
The estimate (17) holds true if a method with optimal computational complexity is used to solve the auxiliary linear systems that appear in (15) or (16). In the multidimensional case, this means that such a preconditioned iterative solver is applied. In the numerical tests presented in our previous papers [9,11,23], the Boomer-AMG implementation from HYPRE [14] of the algebraic multigrid (AMG) preconditioner in the PCG framework is used for this purpose.
The computational complexity depends on the target accuracy of the numerical solution. Assume that the matrix A α corresponds to FEM discretization of the fractional diffusion problem. The rigorous error analysis is a task beyond the standard theory of FEM. Such error bounds are derived in [24] under certain assumptions that are an interplay among data regularity, regularity pickup of the differential operator (e.g., the Laplacian), and the fractional power α . If linear finite elements are used and if f L 2 ( Ω ) , the L 2 error is essentially
| | u u h | | L 2 = O ( h 2 α ) .
The case where the finite difference method (FDM) is applied for discretization of the fractional diffusion problem is analyzed in [10]. The presented error estimates are based on certain equivalence between the lumped mass FEM approximation and the ( 2 d + 1 ) -stencil FDM approximation.
The further analysis integrates the FEM/FDM error estimate (18) in the complexity estimate (17). Assuming the mesh is quasi uniform, N n d , with n h 1 , and therefore
h 2 α = O ( N 2 α / d ) .
Then, for a given mesh parameter h, the degree of BURA approximation k is chosen such that the contributions of the discretization error (18) and the BURA error (14) to the total error are balanced. In this way, and taking into account (17), the following almost optimal computational complexity
N α , k B U R A = O ( N ( log N ) 2 )
is derived in [10].
More recently, the so-called reduced additive BURA-AR and multiplicative BURA-MR methods were introduced in [9]. This study follows the observation that the matrices of part of the auxiliary systems possess very different properties. An important observation is that under certain assumptions, the computational complexity of the reduced BURA can be improved to O ( N log N ) .

4. Bura Preconditioning: Condition Number Estimates

The iterative solution of fractional diffusion problems can be applied in very limited cases. Note that even if a good preconditioner is available, the matrix vector multiplication with A α remains an open problem in the general case.
The study in this paper is motivated by the application of the BURA approximation of A α in the spirit of block-diagonal preconditioners (6), (8) of the coupled problems (4), (7), respectively. In these specific cases, the preconditioned iterative solvers are applied to sparse matrices, and therefore the matrix vector multiplication does not create computational difficulties.
Consider the matrix A α , α ( 0 , 1 ) . Here, we assume the A is a N × N sparse SPD matrix obtained after FDM or FEM discretization of a given elliptical (say Laplacian) boundary value problem with a mesh parameter h. Denote by S p ( A ) = { λ i , h , ψ i , h } i = 1 N the spectrum of A , where { λ i , h } i = 1 N are the monotonically increasing positive eigenvalues and { ψ i , h } i = 1 N are the corresponding normalized eigenvectors. Then
S p ( A α ) = { λ i , h α , ψ i , h } i = 1 N and S p ( A α ) = { λ i , h α , ψ i , h } i = 1 N .

4.1. Bura Preconditioner C α , k B U R A : α ( 0 , 1 )

Following the BURA approximation of A α , we introduce the preconditioner C α , k B U R A of A α by its inverse
C α , k B U R A 1 = λ 1 , h α r α , k ( λ 1 , h A 1 ) .
Lemma 1.
Assume that the degree k is large enough so that E α , k κ α ( A ) < 1 . Then, the following condition number estimate for the BURA preconditioner C α , k B U R A , α ( 0 , 1 ) , holds true
κ C α , k B U R A 1 A α 1 + E α , k κ α ( A ) 1 E α , k κ α ( A ) .
Proof. 
It follows from (19) and (20) that
S p C α , k B U R A 1 A α = λ 1 , h α r α , k ( λ 1 , h / λ i , h ) λ i , h α , ψ i , h i = 1 N ,
and therefore
κ C α , k B U R A 1 A α = max i λ 1 , h α r α , k ( λ 1 , h / λ i , h ) λ i , h α min i λ 1 , h α r α , k ( λ 1 , h / λ i , h ) λ i , h α = max i r α , k ( ζ i , h 1 ) ζ i , h α min i r α , k ( ζ i , h 1 ) ζ i , h α ,
where ζ i , h = λ i , h / λ 1 , h 1 , κ ( A ) for all i = 1 , 2 , , N . Now, applying (13) for ζ i , h 1 [ κ ( A ) 1 , 1 ] [ 0 , 1 ] , and using the equality
r α , k ζ i , h 1 ζ i , h α = 1 + r α , k ζ i , h 1 ζ i , h α ζ i , h α ,
we derive the estimates
1 E α , k ζ i , h α r α , k ζ i , h 1 ζ i , h α 1 + E α , k ζ i , h α , i = 1 , 2 , , N .
Therefore
max i r α , k ( ζ i , h 1 ) ζ i , h α min i r α , k ( ζ i , h 1 ) ζ i , h α max i 1 + E α , k ζ i α min i 1 E α , k ζ i α 1 + E α , k κ α ( A ) 1 E α , k κ α ( A ) ,
which completes the proof. □

4.2. Bura Preconditioner C β , k B U R A : β ( 1 , 0 )

In this case, 1 + β ( 0 , 1 ) . Thus, we consider the preconditioner C β , k B U R A in the form
C β , k B U R A 1 = A λ 1 , h ( 1 + β ) r 1 + β , k ( λ 1 , h A 1 ) .
Lemma 2.
Assume that the degree k is large enough so that E 1 + β , k κ 1 + β ( A ) < 1 . Then, the following condition number estimate for the BURA preconditioner C β , k B U R A , β ( 1 , 0 ) , holds true
κ C β , k B U R A 1 A β 1 + E 1 + β , k κ 1 + β ( A ) 1 E 1 + β , k κ 1 + β ( A ) .
Proof. 
According to (19), the spectrum of the new preconditioner is as follows
S p C β , k B U R A 1 A β = λ 1 , h ( 1 + β ) r 1 + β , k ( λ 1 , h / λ i , h ) λ i , h 1 + β , ψ i , h i = 1 N .
Following the notations and arguments from the proof of Lemma 1, we obtain
κ C β , k B U R A 1 A β = max i r 1 + β , k ( ζ i , h 1 ) ζ i , h 1 + β min i r 1 + β , k ( ζ i , h 1 ) ζ i , h 1 + β .
Hence
κ C β , k B U R A 1 A β max i 1 + E 1 + β , k ζ i , h 1 + β min i 1 E 1 + β , k ζ i , h 1 + β 1 + E 1 + β , k κ 1 + β ( A ) 1 E 1 + β , k κ 1 + β ( A ) ,
and the proof is completed.
It is worth noting that the estimates (21) and (24) are obtained without any regularity assumptions of the related fractional diffusion problem. □

5. Preconditioning of the Coupled Problems: Condition Number Estimates and Computational Complexity

Let us assume that proper FDM or FEM is applied for numerical solution of the coupled problems (4) and (7). Thus, we obtain the block-diagonal preconditioner
C D S = A Δ + D T t T T t K 1 ( A · + I ) I K I ( A Δ + I ) 1 / 2 ,
corresponding to (6) and
C 3 1 = A Δ + I A Δ + I A Δ 0.14 ,
corresponding to (8). As in the continuous case, the preconditioners C D S and C 3 1 provide optimal condition number estimates of the saddle-point matrices corresponding to (4) and (7).
Now, the block-diagonal preconditioners C D S B U R A and C 3 1 B U R A are introduced by replacing the last (interface) blocks in (26) and (27) by their rational approximations, as defined in the previous section. In this way, we obtain
C D S B U R A = A Δ + D T t T T t K 1 ( A · + I ) I K I C 0.5 , k B U R A [ ( A Δ + I ) ] ,
where C 0.5 , k B U R A [ · ] is determined following (20), and
C 3 1 B U R A = A Δ + I A Δ + I C 0.14 , k B U R A [ A Δ ] ,
with C 0.14 , k B U R A [ · ] defined through (23).
The following assumptions will be used in the analysis of the introduced BURA preconditioners:
(A1)
Quasi uniform meshes with mesh parameter h are used for FDM/FEM approximation of the considered coupled problems. Thus, the number of unknowns for the discrete Darcy–Stokes problem is N D S = O ( h 2 ) , and similarly, for the 1D–3D coupling problem, N 3 1 = O ( h 3 ) . In both cases, the number of unknowns related to the interface Γ is N Γ = O ( h 1 ) .
(A2)
Solvers of optimal computational complexity are applied for the systems with sparse SPD matrices which appear in the implementation of the BURA preconditioners.

5.1. Preconditioner C D S B U R A

Lemma 3.
Let the assumptions A1 and A2 be fulfilled, and let k be the minimal degree of the rational approximation of the last block in (28) such that
E 0.5 , k κ 0.5 ( A Δ + I ) c d s ,
where 0 < c d s < 1 is an a priori chosen constant. Then,
κ [ C D S B U R A ] 1 C D S 1 + c d s 1 c d s .
Proof. 
The difference between block-diagonal preconditioners C D S B U R A and C D S is in last blocks only. Consequently, the estimate (30) follows directly from Lemma 1. □
For example, following the estimate from Lemma 3, and choosing c d s = 1 / 2 , we obtain the bound κ [ C D S B U R A ] 1 C D S 3 . Obviously, this condition number can be arbitrary close to 1 if c d s is small enough. However, this might be computationally expensive, taking into account that smaller c d s gives rise to larger degree k of the applied BURA approximation.
The block-diagonal preconditioner C D S B U R A has an optimal convergence rate, provided the assumptions of Lemma 3 are met. This means that O ( 1 ) iterations are sufficient to solve the discrete Darcy–Stokes problem. Therefore, the question is what the computational complexity of solving a system with C D S B U R A is. It can be written in the form
N [ C D S B U R A ] 1 v = N 1 + N 2 + N 3 + N 4 + N 5 ,
where N i is the complexity of solving a system with the i-th diagonal block. It follows from the assumptions (A1) and (A2) that
N 1 + N 2 + N 3 + N 4 = O ( h 2 ) ,
and
N 5 = N [ C 0.5 , k B U R A ( ( A Δ + I ) 0.5 ) ] 1 w = O ( k h 1 ) .
Substituting (13) into (29) and taking into account the sharpness of (13), we obtain
8 e 2 π k / 2 κ ( A Δ + I ) c d s .
Let us recall that κ ( A Δ + I ) = O ( h 2 ) . Thus, the assumptions of Lemma 3 are valid if and only if k = O ( log 2 ( h 1 ) ) , and therefore
N 5 = O ( h 1 log 2 ( h 1 ) ) .
The presented analysis is summarized in the following theorem.
Theorem 1.
Under the assumptions of Lemma 3, the preconditioned iterative solution of a discrete Darcy–Stokes problem with the preconditioner C D S B U R A has an optimal computational complexity.
Proof. 
The statement follows directly from the above estimates of N i , i = 1 , 2 , , 5 , and their sum, that is
N [ C D S B U R A ] 1 v = O ( h 2 ) + O ( h 1 log 2 ( h 1 ) ) =
O ( N D S ) + O ( N D S 1 / 2 log 2 N D S ) = O ( N D S ) .
 □

5.2. Preconditioner C 3 1 B U R A

Lemma 4.
Let the assumptions A1 and A2 be fulfilled. We also assume that k is the minimal degree of the rational approximation such that
E 0.86 , k κ 0.86 ( A Δ ) c 3 1 ,
where 0 < c 3 1 < 1 is a chosen constant. Then,
κ [ C 3 1 B U R A ] 1 C 3 1 1 + c 3 1 1 c 3 1 .
Proof. 
Again, the difference between the block-diagonal preconditioners C 3 1 B U R A and C 3 1 is in the interface blocks only, and the proof is analogous to the one of Lemma 3. □
The block-diagonal preconditioner C 3 1 B U R A has an optimal convergence rate if the assumptions of Lemma 4 are met. This raises the important question about the computational complexity of solving a system with C 3 1 B U R A . In this case, we write it in the form
N [ C 3 1 B U R A ] 1 v = N 1 + N 2 + N 3 .
It follows from the assumptions (A1) and (A2) that
N 1 + N 2 = O ( h 3 ) ,
and
N 3 = N [ C 0.14 , k B U R A ( A Δ 0.14 ) ] 1 w = O ( k h 1 ) .
The analysis of the degree k from Section 5.1 is valid here as well, that is
k = O log 2 ( h 1 ) .
Then, the analogue of Theorem 1 reads as follows.
Theorem 2.
Under the assumptions of Lemma 4, the preconditioned iterative solution of a discrete 3D–1D coupled problem with the preconditioner C 3 1 B U R A has an optimal computational complexity.
Proof. 
The result follows directly from the estimates of N i , i = 1 , 2 , 3 , and their sum as follows
N [ C 3 1 B U R A ] 1 v = O ( h 3 ) + O ( h 1 log 2 ( h 1 ) ) =
O ( N 3 1 ) + O ( N 3 1 1 / 3 log 2 N 3 1 ) = O ( N 3 1 ) .
 □
Theorems 1 and 2 show the asymptotic optimality of the proposed BURA preconditioners. The efficient implementation of the corresponding iterative solution methods depends on the condition number estimates given by Lemmas 3 and 4. They are controlled by the constants c D S and c 3 1 . Some practical aspects related to the balance of these parameters are presented in the next section.
For large discrete linear systems approximating discrete fractional elliptic equations, parallel multigrid solvers are used. The efficiency of these solvers is studied in [25] and the references therein. We refer to recommendations given in this paper.

6. Behavior of the Condition Numbers of the BURA-Based Preconditioners

In this section, we will numerically analyze the behavior of the estimates for the condition numbers of the BURA-based preconditioners. The theoretical estimates derived in Lemmas 1 and 2 will be compared with those that follow from the definitions of the corresponding rational approximations.

6.1. Analysis of κ [ C 0.5 , k B U R A ] 1 A 0.5

From the definition (20) of the preconditioner C α , k B U R A for α ( 0 , 1 ) , we obtain the estimate
κ C 0.5 , k B U R A 1 A 0.5 = max i r 0.5 , k ( ζ i , h 1 ) ζ i , h 0.5 min i r 0.5 , k ( ζ i , h 1 ) ζ i , h 0.5 max η [ 1 , δ ] r 0.5 , k ( η 1 ) η 0.5 min η [ 1 , δ ] r 0.5 , k ( η 1 ) η 0.5 ,
where ζ i , h = λ i , h / λ 1 , h and δ = κ ( A ) . The results of the performed tests are presented in Figure 3.
The data obtained via (32) are denoted as computed estimates, while the theoretical estimates are obtained by applying the upper bound (21) from Lemma 1. For the computed estimates, the BRASIL software [26] has been used. The figure includes theoretical estimates for large enough k for which the condition (29) from Lemma 3 is satisfied.
The first observation is that the computed and theoretical estimates are practically the same, being almost equal to one for k 9 and all considered cases. Of course, it is not necessary to have a relative condition number that is so close to one. Thus, for κ ( A ) = 10 6 , we obtain approximately that κ [ C 0.5 , k B U R A ] 1 A 0.5 3 when k 3 . For κ ( A ) = 10 8 , the same holds true when k 5 .
The relative condition number κ [ C 0.5 , k B U R A ] 1 A 0.5 depends on κ ( A ) . However, the presented results can be assessed as very promising. Assume that for the problem under consideration in Ω R 2 , a quasi-uniform triangulation with mesh parameter h is used. Then, the size of the discrete problem is N = O ( h 2 ) = O ( κ ( A ) ) . Thus, κ ( A ) = 10 8 means N = O ( 10 8 ) , which is more or less the largest size of the most commonly solved 2D problems. In this sense, the tests in Figure 3 are fully representative.

6.2. Analysis of κ [ C 0.14 , k B U R A ] 1 A 0.14

Here, we adopt the approach, assumptions, and notations form the previous subsection. Now, from the definition (23) of the preconditioner C β , k B U R A for β ( 1 , 0 ) , we obtain the estimate
κ C 0.14 , k B U R A 1 A 0.14 = max i r 0.86 , k ( ζ i , h 1 ) ζ i , h 0.86 min i r 0.86 , k ( ζ i , h 1 ) ζ i , h 0.86 max η [ 1 , δ ] r 0.86 , k ( η 1 ) η 0.86 min η [ 1 , δ ] r 0.86 , k ( η 1 ) η 0.86 .
The results shown in Figure 4 are obtained using the estimates (33) and (24).
In this figure, the data for the theoretical estimates start from the first k for which the condition (31) in Lemma 4 is satisfied. Similar to Figure 3, the relative condition number is close to one if k 10 for both values of κ ( A ) { 10 6 , 10 8 } . Here, for κ ( A ) = 10 6 , we see that κ [ C 0.14 , k B U R A ] 1 A 0.14 3 when k 5 . For κ ( A ) = 10 8 , the same holds true when k 8 .
Assume now that for the problem of coupling a 1D structure embedded in Ω R 3 , a quasi-uniform mesh with a parameter h is used. Then, the size of the discrete problem is N = O ( h 3 ) and therefore κ ( A ) = 10 8 means N = O ( 10 12 ) . Such extremely large-scale 3D problems are a challenge for even the most powerful supercomputers today. Thus, κ ( A ) = 10 8 can be considered as a limit case when analyzing the preconditioner C 0.14 , k B U R A .

6.3. Comparative Summary

The preconditioned interface problems discussed in the previous two sections are significantly different, due to the alternating sign of the degrees α and β . However, the condition number estimates for the proposed BURA preconditioners are very similar. Thus, the definition (23) implies the equality
κ C 0.14 , k B U R A 1 A 0.14 = κ C 0.86 , k B U R A 1 A 0.86 ,
therefore the analysis can be reduced to the case of positive fractional powers. From this point of view, the estimates (21) from Lemma 1 and (24) from Lemma 2 show certain advantages of C 0.86 , k B U R A compared to C 0.5 , k B U R A . This follows from the better exponential accuracy E α , k of the BURA method for larger α ( 0 , 1 ) . However, this conclusion is asymptotically true, that is, for large enough k. The data in Figure 3 and Figure 4 show a different picture. This is the motivation to present more detailed information about the behavior of κ C 0.5 , k B U R A 1 A 0.5 and κ C 0.14 , k B U R A 1 A 0.14 for 3 k 9 in Table 1.
The presented data confirm the promising assessment of the BURA preconditioners. We observe a stable behavior of the computed estimates of the condition numbers. If needed, the tabular results can be used for good enough piecewise linear interpolation for arbitrary δ ( 10 5 , 10 8 ) .
As noted in the Introduction, our study is partly motivated by the multigrid methods for discrete fractional Sobolev spaces presented in [5]. For the 1D fractional Laplacian in the interval [ 0 , 1 ] and uniform mesh, the reported relative condition number of the multigrid preconditioner for α = 0.5 and for the maximal size n = 512 is κ 0.5 M G = 3.2 . The spectrum of A for this model problem is known. Thus, to compare κ 0.5 M G = 3.2 with the data from Table 1, we recall that n = 512 corresponds to κ ( A ) 10 5 . In the case of negative fractional power, we interpolate the presented tabular data to obtain the corresponding relative condition number κ 0.14 M G 16 . These simple examples illustrate the competitiveness of the BURA preconditioners. The potential benefits of the new method will be even greater for coupled problems with more complex geometry Γ .

7. Concluding Remarks

The fractional diffusion operators are non-local, and the related matrices are dense. Furthermore, when using the spectral definition, the matrix A α is generally non-computable. In such a situation, we are not even able to compute the residual, which should directly rule out the idea of applying any iterative solution methods to systems in the form A α u = f .
Why, then, do we discuss preconditioning methods for A α ? For the considered coupled problems such preconditioning makes sense, as the corresponding systems have sparse matrices (so computing of the residue is not a problem), while the fractional power matrices appear only in the block-diagonal preconditioners (6) and (8).
Our study is motivated by the results of the recently published papers [5,13,17], where the developed multigrid methods for discrete fractional Sobolev spaces are applied to the Darcy–Stokes and 3D–1D coupled problems. Alternatively, we have proposed a fully algebraic preconditioning approach.
For α ( 0 , 1 ) , the new preconditioner C α , k B U R A relies directly on the best uniform rational approximation (BURA) of degree k of A α . An important further contribution is the construction of C β , k B U R A for β ( 1 , 0 ) leading to the equality
κ C β , k B U R A 1 A β = κ C 1 + β , k B U R A 1 A 1 + β .
This allowed us to develop a unified approach to the relative condition number estimates (21) and (24), applied later in Theorem 1 and Theorem 2, where the optimal computational complexity for both considered coupled problems is proved. The numerical study of κ [ C 0.5 , k B U R A ] 1 A 0.5 and κ [ C 0.14 , k B U R A ] 1 A 0.14 is also important for the value of the paper. As shown, a very small degree of k is practically sufficient to ensure high quality BURA preconditioning.
In general, it is clear that a required accuracy can be obtained by applying different known methods for spectral fractional diffusion problems. Here, we apply the BURA method, motivated by its advantages, as proven in [19] (see the discussion at the beginning of Section 3).
For the simplest case of Γ = [ 0 , 1 ] , the comparison with the multigrid preconditioners from [5] illustrates the competitiveness of the newly proposed C 0.5 , k B U R A and C 0.14 , k B U R A . For most of the BURA preconditioners, even stronger numerical advantages can be expected when more complex coupled problems are considered, due to the more general scope of applicability of the presented results. Among the other advantages, we note the following: BURA preconditioners are purely algebraic, so regularity assumptions are not made in the theoretical analysis. The theoretical estimates (21) and (24), as well as those shown in Table 1, are exactly the same for every type of geometry of Γ . The BURA approach and the obtained condition number estimates are applicable to a wider class of problems, including, for example, coupling of a 2D manifold embedded in Ω R 3 .

Author Contributions

Conceptualization, S.M.; Formal analysis, I.L.; Funding acquisition, S.H. and S.M.; Methodology, S.H. and S.M.; Project administration, S.M.; Software, I.L.; Visualization, I.L.; Writing—original draft, S.M.; Writing—review & editing, S.H. and I.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Grant No. BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014–2020) and co-financed by the European Union through the European Structural and Investment and by the Bulgarian National Science Fund under grant No. DFNI-DN12/1.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We acknowledge the provided access to the e-infrastructure and support of the Centre for Advanced Computing and Data Processing, with the financial support by the Grant No. BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014–2020) and co-financed by the European Union through the European Structural and Investment funds. The work is partially supported by the Bulgarian National Science Fund under grant No. DFNI-DN12/1.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Keyes, D.; McInnes, L.; Woodward, C.; Gropp, W.; Myra, E.; Pernice, M.; Bell, J.; Brown, J.; Clo, A.; Connors, J.; et al. Multiphysics simulations: Challenges and opportunities. Int. J. High Perform. Comput. Appl. 2013, 27, 4–83. [Google Scholar] [CrossRef] [Green Version]
  2. Mehdizadeh Khalsaraei, M.; Shokri, A.; Molayi, M. The new high approximation of stiff systems of first order IVPs arising from chemical reactions by k-step L-stable hybrid methods. Iran. J. Math. Chem. 2019, 10, 181–193. [Google Scholar]
  3. Nordsletten, D.; Niederer, S.; Nash, M.; Hunter, P.; Smith, N. Coupling multi-physics models to cardiac mechanics. Prog. Biophys. Mol. Biol. 2011, 104, 77–88. [Google Scholar] [CrossRef]
  4. Shokri, A.; Mehdizadeh Khalsaraei, M.; Atashyar, A. A new two-step hybrid singularly P-stable method for the numerical solution of second-order IVPs with oscillating solutions. Iran. J. Math. Chem. 2020, 11, 113–132. [Google Scholar]
  5. Bærland, T.; Kuchta, M.; Mardal, K.A. Multigrid methods for discrete fractional Sobolev spaces. Siam J. Sci. Comput. 2019, 41, A948–A972. [Google Scholar] [CrossRef]
  6. Kuchta, M.; Mardal, K.A.; Mortensen, M. Preconditioning trace coupled 3D-1D systems using fractional Laplacian. Numer. Methods Partial. Differ. Equ. 2019, 35, 375–393. [Google Scholar] [CrossRef] [Green Version]
  7. Layton, W.J.; Schieweck, F.; Yotov, I. Coupling fluid flow with porous media flow. Siam J. Numer. Anal. 2002, 40, 2195–2218. [Google Scholar] [CrossRef]
  8. 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]
  9. 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]
  10. 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] [Green Version]
  11. 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]
  12. Galvis, J.; Sarkis, M. Non-matching mortar discretization analysis for the coupling Stokes-Darcy equations. Electron. Trans. Numer. Anal 2007, 26, 350–384. [Google Scholar]
  13. Holter, K.E.; Kuchta, M.; Mardal, K.A. Robust preconditioning for coupled Stokes–Darcy problems with the Darcy problem in primal form. Comput. Math. Appl. 2021, 91, 53–66. [Google Scholar] [CrossRef]
  14. Falgout, R.D.; Yang, U.M. hypre: A Library of High Performance Preconditioners. In Computational Science—ICCS 2002; Sloot, P.M.A., Hoekstra, A.G., Tan, C.J.K., Dongarra, J.J., Eds.; Lecture Notes in Computer Science; Springer: Berlin/Heidelberg, Germany, 2002; Volume 2331, pp. 632–641. [Google Scholar]
  15. Zhang, H.; Guo, J.; Lu, J.; Li, F.; Xu, Y.; Downar, T. An assessment of coupling algorithms in HTR simulator TINTE. Nucl. Sci. Eng. 2018, 190, 287–309. [Google Scholar] [CrossRef]
  16. Holter, K.E.; Kuchta, M.; Mardal, K.A. Robust Preconditioning of Monolithically Coupled Multiphysics Problems. 2020. Available online: https://arxiv.org/abs/2001.05527 (accessed on 28 February 2020).
  17. Kuchta, M.; Laurino, F.; Mardal, K.A.; Zunino, P. Analysis and approximation of mixed-dimensional PDEs on 3D-1D domains coupled with Lagrange multipliers. Siam J. Numer. Anal. 2021, 59, 558–582. [Google Scholar] [CrossRef]
  18. Duan, B.; Lazarov, R.; Pasciak, J. Numerical approximation of fractional powers of elliptic operators. IMA J. Numer. Anal. 2020, 40, 1746–1771. [Google Scholar] [CrossRef] [Green Version]
  19. Hofreither, C. A Unified View of Some Numerical Methods for Fractional Diffusion. Comput. Math. Appl. 2020, 80, 332–350. [Google Scholar] [CrossRef]
  20. Stahl, H. Best uniform rational approximation of xα on [0, 1]. Bull. Am. Math. Soc. 1993, 28, 116–122. [Google Scholar] [CrossRef] [Green Version]
  21. Saff, E.B.; Stahl, H. Asymptotic distribution of poles and zeros of best rational approximants to xα on [0, 1]. In Topics in Complex Analysis; Banach Center Publications; Institute of Mathematics, Polish Academy of Sciences: Warsaw, Poland, 1995; Volume 31. [Google Scholar]
  22. Hofreither, C. An algorithm for best rational approximation based on barycentric rational interpolation. Numer. Algorithms 2021, 88, 365–388. [Google Scholar] [CrossRef]
  23. Harizanov, S.; Kosturski, N.; Margenov, S.; Vutov, Y. Neumann fractional diffusion problems: BURA solution methods and algorithms. Math. Comput. Simul. 2021, 189, 85–98. [Google Scholar] [CrossRef]
  24. Bonito, A.; Pasciak, J. Numerical approximation of fractional powers of elliptic operators. Math. Comput. 2015, 84, 2083–2110. [Google Scholar] [CrossRef] [Green Version]
  25. Čiegis, R.; Starikovičius, V.; Margenov, S.; Kriauzienė, R. Scalability analysis of different parallel solvers for 3D fractional power diffusion problems. Concurr. Comput. Pract. Exp. 2019, 31, e5163. [Google Scholar] [CrossRef]
  26. Software BRASIL. Available online: https://baryrat.readthedocs.io/en/latest/#baryrat.brasil (accessed on 28 February 2020).
Figure 1. Darcy–Stokes system in a rectangle domain Ω = Ω S Ω D : model problem with simple geometry (left) and more realistic problem with a more general geometry of the interface (right).
Figure 1. Darcy–Stokes system in a rectangle domain Ω = Ω S Ω D : model problem with simple geometry (left) and more realistic problem with a more general geometry of the interface (right).
Mathematics 10 00780 g001
Figure 2. Coupled 3D–1D systems in a hexahedral domain Ω : model problems with single line geometry of the interfaces (left) and more realistic problem with more general geometry of the interface (right).
Figure 2. Coupled 3D–1D systems in a hexahedral domain Ω : model problems with single line geometry of the interfaces (left) and more realistic problem with more general geometry of the interface (right).
Mathematics 10 00780 g002
Figure 3. Theoretical estimates (21) vs. computed estimates (22) of κ [ C 0.5 , k B U R A ] 1 A 0.5 with respect to the BURA degree k, for δ = κ ( A ) { 10 6 , 10 8 } .
Figure 3. Theoretical estimates (21) vs. computed estimates (22) of κ [ C 0.5 , k B U R A ] 1 A 0.5 with respect to the BURA degree k, for δ = κ ( A ) { 10 6 , 10 8 } .
Mathematics 10 00780 g003
Figure 4. Theoretical estimates (24) vs. computed estimates (25) of κ [ C 0.14 , k B U R A ] 1 A 0.14 with respect to the BURA degree k, for δ = κ ( A ) { 10 6 , 10 8 } .
Figure 4. Theoretical estimates (24) vs. computed estimates (25) of κ [ C 0.14 , k B U R A ] 1 A 0.14 with respect to the BURA degree k, for δ = κ ( A ) { 10 6 , 10 8 } .
Mathematics 10 00780 g004
Table 1. Computed estimates of κ [ C 0.5 , k B U R A ] 1 A 0.5 and κ [ C 0.14 , k B U R A ] 1 A 0.14 with respect to 3 k 9 , for δ = κ ( A ) { 10 5 , 10 6 , 10 7 , 10 8 } .
Table 1. Computed estimates of κ [ C 0.5 , k B U R A ] 1 A 0.5 and κ [ C 0.14 , k B U R A ] 1 A 0.14 with respect to 3 k 9 , for δ = κ ( A ) { 10 5 , 10 6 , 10 7 , 10 8 } .
k α = 0.5 β = 0.14
δ = 10 5 δ = 10 6 δ = 10 7 δ = 10 8 δ = 10 5 δ = 10 6 δ = 10 7 δ = 10 8
31.463.239.9831.503.7423.69168.621221.51
41.381.463.2910.191.385.6437.57270.26
51.081.431.463.781.091.9610.2071.12
61.031.081.461.661.061.133.4321.38
71.021.061.181.461.011.091.587.31
81.011.031.081.341.011.051.092.95
91.001.011.031.081.001.011.091.54
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Harizanov, S.; Lirkov, I.; Margenov, S. Rational Approximations in Robust Preconditioning of Multiphysics Problems. Mathematics 2022, 10, 780. https://doi.org/10.3390/math10050780

AMA Style

Harizanov S, Lirkov I, Margenov S. Rational Approximations in Robust Preconditioning of Multiphysics Problems. Mathematics. 2022; 10(5):780. https://doi.org/10.3390/math10050780

Chicago/Turabian Style

Harizanov, Stanislav, Ivan Lirkov, and Svetozar Margenov. 2022. "Rational Approximations in Robust Preconditioning of Multiphysics Problems" Mathematics 10, no. 5: 780. https://doi.org/10.3390/math10050780

APA Style

Harizanov, S., Lirkov, I., & Margenov, S. (2022). Rational Approximations in Robust Preconditioning of Multiphysics Problems. Mathematics, 10(5), 780. https://doi.org/10.3390/math10050780

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