Matrices associated to two conservative discretizations of Riesz fractional operators and related multigrid solvers

In this article, we focus on a two‐dimensional conservative steady‐state Riesz fractional diffusion problem. As is typical for problems in conservative form, we adopt a finite volume (FV)‐based discretization approach. Precisely, we use both classical FVs and the so‐called finite volume elements (FVEs). While FVEs have already been applied in the context of fractional diffusion equations, classical FVs have only been applied in first‐order discretizations. By exploiting the Toeplitz‐like structure of the resulting coefficient matrices, we perform a qualitative study of their spectrum and conditioning through their symbol, leading to the design of a second‐order FV discretization. This same information is leveraged to discuss parameter‐free symbol‐based multigrid methods for both discretizations. Tests on the approximation error and the performances of the considered solvers are given as well.

Here, we focus on two different finite volume-based discretizations of the following two-dimensional boundary-value steady-state conservative Riesz fractional diffusion equation (FDE) of order 2 − , 2 − , with 0 < , < 1 and with absorbing boundary conditions: where The outline of this article is the following. In Section 2, we introduce some preliminary tools. Section 3 contains the FVE discretization of (1). In Section 4, we provide the FV discretization of (1) and the spectral analysis of the corresponding coefficient matrix, which is our main contribution. In Section 5, we define an ad hoc multigrid solver and a banded preconditioner, and in Section 6, we report some numerical results. Our conclusions are drawn in Section 7.

PRELIMINARIES
This section contains various preliminaries on fractional derivatives (Section 2.1) and Toeplitz matrices (Section 2.2) needed in the rest of this article.

Fractional operators
The Riesz fractional operator in the x-variable is defined as where the left (+) and right (−) derivatives are given in the RL form, that is, with Γ(⋅) being the gamma function. Similarly one can define the Riesz fractional operator in the y-variable.
An alternative definition of the left and right fractional derivatives is based on the shifted Grünwald-Letnikov (GL) formulas. For any p, q ∈ Z and 0 < h x ≪ 1, we define the p-shifted and q-shifted GL fractional derivatives of order 1 − , with 0 < < 1, as 1− u(x, y) GL,p,h x + 1− u(x, y) where t (1− ) . Again similar definitions can be given for the fractional derivatives in the y-variable.
Remark 1. Consider u ∶ R 2 → R with supp(u) ∈ [0, 1] 2 and the equispaced grid with N x , N y ∈ N. Then Equation (2) with p = 0 can be written as 1− u(x i , y j ) GL,0,h x + Matrix G +,0 is a Toeplitz matrix and represents the left fractional derivative operator. The choice of p ≠ 0 yields the same structured matrix, but with the diagonals shifted to the right of p positions, if p > 0, and the diagonals shifted to the left of |p| positions, if p < 0. We denote such an operator by G +,p . Note that in the case of p > 0 we have to compute p new coefficients t (1 (1− ) N x +p−1 to fill the bottom left diagonals. Furthermore, when q = 0, the right fractional derivative operator in Equation (3) is G −,0 = −G T +,0 . If q ≠ 0 then we denote such an operator by G −,q and it holds G −,q = −G T +,q . Let h x > 0, then, under proper hypothesis (see Reference 15) it can be proven that, for p 1 , p 2 , q 1 , q 2 ∈ Z with p 1 ≠ p 2 and q 1 ≠ q 2 , 1− u(x, y) RL where w p = 1− −2p 2 2(p 1 −p 2 ) and w q = 1− −2q 2 2(q 1 −q 2 ) with p = (p 1 , p 2 ) and q = (q 1 , q 2 ). Of course, the same reasoning applies to the y-variable as well.
where h x is the step length. Due to the FV discretization (see Section 4), we are required to evaluate the fractional derivative operators between two grid points, for example, in , which leads to This motivates the need of a non-integer shift.
In References 10 and 16, a non-integer shift was used to define a first-order FV approximation. A second-order FV scheme is still missing in the literature. In order to fill this gap in Section 4, we note that the validity of Equation (5) extends also to the case where p 1 , p 2 , q 1 , q 2 ∈ R with p 1 ≠ p 2 and q 1 ≠ q 2 . The reader can easily verify this assuming that p 1 , p 2 , q 1 , q 2 are real and following the same argument of the proof of Theorem 1 in Reference 15.
As a consequence, for a generic step length h x > 0 and for p 1 ≠ p 2 ∈ Z + 1 2 , q 1 ≠ q 2 ∈ Z + 1 2 , we can write We refer the reader to Theorem 1 for the matrix form of Equation (7).

Toeplitz matrices and their symbol
As already observed in Remark 1, when considering uniform meshes, the discretization of the fractional operators yields a Toeplitz matrix. In order to explore the properties of such matrices, we recall some basic definitions, see, for example, Reference 17.

Definition 1.
A Toeplitz matrix T N ∈ C N×N has constant coefficients along the diagonals, namely Remark 3. If T N is generated by f , then T H N is generated by f . Hence, T N is Hermitian whenever f is a real function. Furthermore, in the case where coefficients t k are real it holds that f (x) = f (−x), hence if f is real then it is also even.
In the case of a two-dimensional operator, the discretization of the considered fractional derivative operators yields a two-level Toeplitz matrix, which is a block-Toeplitz with Toeplitz blocks (BTTB) matrix. The BTTB matrix of order N = N x N y generated by f is n i ∈ R n i ×n i has entry (r, s)th equals 1 if s − t = j i and 0 elsewhere. Remark 4. An interesting computational property of the unilevel Toeplitz matrix T N is that the matrix-vector product can be performed in O(N log N) through the fast Fourier transform (FFT) algorithm. 18 A similar property holds for BTTB as well.
Remark 5. The symbol has many properties, one of which is that it allows to estimate the range of the eigenvalues of an Hermitian Toeplitz matrix.
In the following, we aim at writing explicitly the symbol of the (properly scaled) Toeplitz matrices representing the discretized operators in Equation (6) of Remark 2. Such symbol will be useful in the computation of the symbol for the FV discretization of (1) performed in Section 4. Having this in mind, we start with a couple of intermediate results.

Proposition 1.
Let N x ∈ N, then it holds that G +,0 , G −,0 defined in Equation (4) are such that Proof. According to the definition of symbol and by means of the generalized Newton binomial, it holds which completes the proof. ▪ The following result, proved in Appendix A, is needed for later analysis.

Lemma 1.
For all x ∈ [0, ] it holds that and As a corollary to Proposition 1, we obtain the symbol of the Toeplitz matrices representing the second-order discretization of the fractional derivative operators given in Equation (5).
, and p, q ∈ Z 2 , then Equation (5) can be written as follows where u (j) is defined in Equation (4) and Proof. According to the definition of symbol, shifting the diagonals by p positions to the right or left consists in multiplying the symbol by e −ipx or e ipx , respectively. Therefore, the proof follows by Remark 1. ▪ We are now ready to provide the symbol of the Toeplitz matrices corresponding to the fractional left and right operators evaluated at the midpoint x i− 1 2 given in Equation (6).
, and p, q ∈ Z 2 + 1 2 , then Equation (7) can be written as follows where u (j) is defined in Equation (4) and H +,p = T N x ) .
Proof. We only provide the proof in the case of the left fractional derivative since for the other case the proof follows the same steps. From Equation (7), which represents the ith row of the matrix-vector product in Equation (10), we have that the resulting matrix is a Toeplitz generated by which completes the proof. ▪

FV-TYPE DISCRETIZATIONS
Here we review the idea of an FV-based discretization applied to problem (1). First, we cover the domain Ω with a mesh ∪ n i=1 Q i , where (Q i ∩ Q j ) = 0, i ≠ j, with the Lebesgue measure. Then, by integrating over the control volumes Q i , the approximation of the differential operator is reduced to the approximation of the fractional derivative of order 1 − , ∈ { , }, on cell boundaries. In particular, in Section 3.1, we recall the FVE approach, where the solution is approximated in the space of C 0 finite elements and then fractionally derived using exact formulas for fractional derivatives of a polynomial. The classical FV approach will be treated in Section 4.
In our specific case, given N x , N y ∈ N, we partition the domain Ω = [a 1 , b 1 ] × [a 2 , b 2 ] into N x × N y equal rectangles. More specifically, letting . By integrating Equation (1) over Q ij , we obtain S 1 + S 2 = S 3 , where We approximate S 3 by means of the tensor product of Simpson's rules, which is an order 3 scheme, so that the approximation of the right-hand side will not influence the solution and the comparison of the FV and FVE discretization approaches. Therefore, Performing the integration in dx, one can rewrite S 1 as At this point, we can proceed in two different ways: 1. Either approximating again S 1 as and u(x, y) by a piecewise polynomial and finally computing exactly the remaining integrals thanks to the exact formulas for the fractional derivatives of polynomials (FVE approach); 2. Or approximating the integrals with the midpoint rule, leading to and using the Grünwald formulas for the point values of the fractional derivatives (FV approach).
The order of accuracy of (13) can be understood by reinterpreting the product h y appearing in (14) as the result of applying the midpoint rule to the integral appearing in (13). Finally, the order of accuracy of the fully discrete scheme obtained is further capped by the choice of polynomial spaces in the FVE approach or the order of the Grünwald formulas in the FV approach.
A similar reasoning of course applies to S 2 .

FVE discretization matrices and their spectral study
FVE have already been applied to FDE problems in References 4,9,and 11. In this approach, the unknown function is sought into a classical tensor product Q1 finite element space on the dual grid whose vertices are at the centers of the control volumes Q ij . In this view, we consider the basis functions for k = 1, … , N x , and define similarly y l (y) with y l in place of x k and h y in place of h x . Then, considering as unknowns the function values u ij at the center of Q ij , we replace u(x, y) in Equation (13) Since the support of y l (y) is compact, then Therefore, the FVE discretization of the FDE problem in (1) yields the linear system where the right-hand side b follows from Equation (11), the solution is u = {u kl } N x ,N y k,l=1 and the coefficient matrix is the following N × N, with N = N x N y , matrix with ) , ) .
The grid dependent scale factors are r = As already observed in Reference 14, in the one-dimensional case with constant diffusion coefficients, the symbol of the coefficient matrix is , which is a nonnegative function with a unique zero of order lower than 2 at x = 0. In the case of a two-dimensional equation with constant diffusion coefficients the symbol of A FVE iŝ Remark 6. Note thatĝ 2D (x, y) has a unique zero of order lower than 2 at (x, y) = (0, 0). This is because the symbol of the mass matrix is a strictly positive function.

FV DISCRETIZATION MATRICES AND THEIR SPECTRAL STUDY
FV discretizations consider as unknowns the point values of the function u ij at the centers of the control volumes. Differently from the FVE approach, after having integrated, the fractional derivatives of order 1 − , ∈ { , }, on the boundary of each control volume are now approximated by a fractionally shifted Grünwald formula; by choosing half integer shifts, these fractional derivatives are expressed in terms of the unknowns u ij . First-order accurate FV discretizations for FDE problems appeared in References 5, 10, and 16. Here we build a second-order scheme by imposing some reasonable constraints on the shift parameters involved in the approximation of the fractional derivatives. In addition, on the same line of what has been done in Reference 14, we provide a spectral study of the resulting coefficient matrices which allows to build ad-hoc solvers for the associated linear systems in Section 5.
Let us go back to (13). The choice of approximating S 1 as in (14) and using Equation (10) yields a N × N linear system, whose structure of the coefficient matrix A FV , except for the mass matrices that are replaced by identities, is the same as A FVE in Equation (16). In detail, we have to solve where b follows from Equation (11) while M ,R is obtained by M ,L shifting its diagonals one position forward, that is, . The matrices M ,L , M ,R are similarly defined.

Properties of the symbol of A FV
In the following, we study the properties of A FV and we explain what is a good choice for the shifting parameters with Having in mind the design of an ad-hoc multigrid method for the linear systems associated to A FV , we ask that F (p,q) (x) is a nonnegative function with a unique zero (see Section 5.1 for more details). Let us first require that F (p,q) (x) is a real-valued function. Since there are many free parameters we fix q = p. Under this constraint function F (p,p) (x) reads as which is zero ∀x ∈ (− , ] and ∀p 1 , p 2 , and this implies that F (p,p) (x) is a real-valued function independently of p.
In order to make a reasonable choice of p, we numerically check how the relative 2-norm approximation error varies with p while solving (17) in the case where K x = K y = 1 and solution u(x, y) with related forcing term f (x, y) are the ones reported in Section 6. Many tests show that choosing p 1 , p 2 too far from 0 leads to an increase in the error. Hence, we fix Figure 1 shows the relative 2-norm error for and varying , . We note that the optimal p seems to be p =  due to the symmetry of formulas in Equation (10) with respect to the shifting parameters p 1 , p 2 and q 1 , q 2 , respectively. Therefore, from now onwards, we will fix p =  The numerical results in Section 6 show that such a choice of p and q leads to a second-order accurate numerical scheme for Equation (1) (see Figure 2). has almost one-third of the approximation error than p = (  . Moreover, some preliminary numerical checks, which are not reported here, seem to indicate that the resulting coefficient matrix is positive definite and therefore it could be another interesting combination to investigate. We now check whether for p = , the symbol F (p,p) (x) is nonnegative with a unique zero. For the sake of readability, we omit the superscript (p, p) in the symbol and rewrite it as Then, from Equation (5), we have w p = 2− 2 and from Lemma 1 and the Euler formulas, we have ) .
The following theorem answers positively to our request of having a symbol M ,L − M ,R which is nonnegative with a single zero. The proof follows from the study of the two multiplicative factors of the symbol and is reported in Appendix B. , which is in accordance with the limit case = 0 where we have F 0 (x) = 2(2 − 2 cos x), that is, a multiple of the Laplacian symbol.
It is easy to see that the properties of F (x) transfer to the symbol of A FV . First recall that F (x) is the symbol of M ,L − M ,R in Equation (19). Therefore, multiplying by the scaling parameters and diffusion coefficients we have Similarly, along the second spatial dimension, A y = sK y T N y (F (y)).

Therefore,
If we suppose r s → c, with c ∈ R + , when N x , N y → ∞, then from Theorem 2 the following corollary immediately follows. Corollary 2. Let , ∈ (0, 1), r s → c as N x , N y → ∞ and take constant diffusion coefficients, then the symbol  , (x) is a nonnegative function that has a unique zero at (x, y) = (0, 0) of order min{2 − , 2 − }.

SYMBOL-BASED FAST SOLVERS
Based on the analysis performed in Section 4, in this section, we propose two iterative strategies for solving (15) and (17). Precisely, we present a multigrid method with damped Jacobi as smoother and a band preconditioner whose inverse is approximated through one iteration of the aforementioned multigrid.

Multigrid methods
Multigrid methods combine two iterative methods known as smoother and coarse grid correction (CGC); for more details see, for example, References 20 and 21. The smoother is typically a simple stationary iterative method. The multigrid algorithm can be figured out starting from the two-grid case. One step of a two-grid method is obtained by: (1) computing an initial approximation by few iterations of a pre-smoother, (2) projecting and solving the error equation into a coarser grid, (3) interpolating the solution of the coarser problem, (4) updating the initial approximation, and finally (5) applying a few iterations of a post-smoother to further improve the approximation. Since the coarser grid could be too large for a direct computation of the solution, the same idea can be recursively applied obtaining the so-called V-cycle method.
A common approach to define the coarser operator, known as geometric approach, consists in rediscretizing the same problem on the coarser grid. This approach has the advantage of maintaining the same structure of the coefficient matrix at each level, allowing fast matrix-vector products exploiting the Toeplitz structure (see Remark 4). On the other hand, the coarser problems need to be properly scaled and the result is usually less robust than the so-called Galerkin approach. The latter, for a given linear system A N x = b, A N ∈ C N×N , defines the coarser matrix as A K = P T N A N P N , where P N ∈ C N×K is the full-rank prolongation matrix, while P T N is the restriction operator. The Galerkin approach is useful for the convergence analysis, but in practice it could be computationally too expensive for FDE problems.
The convergence of the V-cycle relies on the so-called smoothing property and approximation property (see Reference 22). In order to discuss the convergence analysis of V-cycle applied either to (15) or (17), we consider constant diffusion coefficients and weighted Jacobi as smoother. Under these assumptions and because of the Toeplitz structure of the considered matrices, the weighted Jacobi coincides with the relaxed Richardson iteration which is well known to satisfy the smoothing property for positive definite matrices, whenever it is convergent. 23 Moreover, thanks to Remarks 6 and 8 and Theorem 2, the approximation property holds with the same projectors as in the case of the Laplacian (see Reference 24) for both FV and FVE approaches.
Since matrices A FV and A FVE are both sums and products between diagonal matrices and dense block Toeplitz matrices, thanks to Remark 4 the matrix-vector product can be performed in O (N log N) operations, without assembling the coefficient matrix, and the storage only requires O(N) elements.
We stress that, by using a Galerkin approach, the block Toeplitz-like structure of A FV and A FVE at the coarser levels is lost, while implementing the geometric approach allows to perform the matrix-vector products by FFT at each coarser grid. Therefore, our multigrid hierarchy is built through the geometric approach and the amount of levels is given by lvl = ⌊log 2 (N x )⌋, that is, the coarsest level has size 1 × 1. Note that in order to make the V-cycle properly working, the linear systems must be scaled such that the right-hand side does not contain any grid dependent scaling factor. Therefore, we scale both A FV x = b and A FVE x = b by h x h y . Similar scalings of course apply also to all coarser levels.
At each iteration of V-cycle one iteration of relaxed Jacobi as pre-and post-smoother is performed. The relaxation parameter is estimated through the approach introduced in Reference 25. Such estimation is obtained by: (1) rediscretizing Equation (1) over a coarser grid (Ñ x , Ñ y ≤ 2 4 ) and keeping the same scaling factors r, s as in the original coefficient matrix, (2) computing the spectrum of the Jacobi iteration matrix, (3) choosing the weight in such a way that the whole spectrum is contained inside a complex set which is the sum of a semicircle and a line, and is motivated by the need of clustering the spectrum of the Jacobi iteration matrix inside the unitary circle. Note that õ(x) yields a set O that is slightly smaller than the unitary circle in such a way that possible outliers are still smaller than 1 in modulus. Our numerical tests in Section 6 confirm that choosing = 0.4, as done in Reference 25, leads to a linearly convergent algorithm.

Banded preconditioner
In References 2 and 14, it was, respectively, proven that coefficients t (  . This motivates the choice, of a band truncation of the discretized fractional operators. Here we consider a band truncation of matrices G ,N x , G ,N y and M ,L , M ,R for FVE and FV, respectively. The resulting block-banded-banded-block matrix Ã is used as GMRES preconditioner. Instead of inverting Ã, we apply one iteration of V-cycle before each iteration of GMRES. The resulting GMRES preconditioner is denoted by  VB , where B is an odd integer number which denotes the block bandwidth and the bandwidth of each block. We expect that for , ≈ 1, preconditioner  VB will perform better for FV than FVE, due to the almost quadratic decay of the coefficient matrix entries in the FV approach compared to the almost linear one in the FVE case. The hierarchy of  VB is built through the geometric approach. On the other hand, by projecting Ã at the coarser levels, the band structure is preserved and the bandwidth does not grow. Therefore, a more robust approach could be obtained building the hierarchy of  VB by means of the Galerkin approach. However, the loss of the diagonal-times-Toeplitz structure would make it harder to estimate the relaxation parameter of Jacobi and a different smoother should be adopted. Note that, due to the band structure of Ã, the matrix-vector product at each level has a linear cost with respect to the matrix size and, as a consequence, each preconditioning iteration costs O(N).

NUMERICAL RESULTS
In this section, we check the second-order convergence of the FV scheme proposed in Section 4 and we test the performances of the methods presented in Section 5 when applied to both (15) and (17). Precisely, we compare the V-cycle algorithm given in Section 5.1 as both main solver (denoted by V) and GMRES preconditioner (denoted by  V ), with the banded preconditioner  VB given in Section 5. Let us consider functionũ(x) = x 2 (1 − x) 2 , x ∈ Ω = [0, 1]. From Reference 13, the exact Riesz fractional derivatives of order 1 − and 2 − ofũ, are where (a 1 , a 2 , a 3 ) = (2, −12, 24). In the following examples, we consider u(x, y) =ũ(x)ũ(y) and build the exact forcing term v(x, y) through the formulas in Equation (20), for these two choices of the diffusion coefficients: • Choice 1: K x (x, y) = K y (x, y) = 1; • Choice 2: 13 K x (x, y) = K y (x, y) = e 4x+4y .

Example 1.
First, we test the accuracy provided by both FVE and FV approaches while considering Choice 1. Figure 2b reports the relative 2-norm error in FV (E FV ) and in FVE (E FVE ), while Figure 2a reports the ratio between the two as N x increases and , vary.
In Figure 2a, a comparison with the black line representing the square of the step length h x (= h y ) confirms the convergence of order 2 for both FV and FVE.
When the ratio between the errors in Figure 2b is smaller than 1, then the FV approach allows better approximation of the solution than the FVE approach. We note that FV has a lower approximation error than FVE in the cases where , ≤ 0.5. Especially, when , ≈ 1 the error in FVE is decreasing faster than in FV, therefore we expect FVE to yield better results than FV when N x > 2 11 − 1. On the contrary, when , ≈ 0 the error in FV decreases faster and reaches almost half the error of FVE for N x = 2 11 − 1. Therefore, it is reasonable to expect further improvements in approximation error for FV with respect to FVE when N x > 2 11 − 1. Further tests, which are not reported here, show that similar results are achieved also for Choice 2.

Example 2.
We now test the behavior of our proposals for solving the two linear systems obtained from FVE and FV when considering Choice 2. Note that, when considering the simpler case of Choice 1, the overall iterations and CPU times would decrease and the Conjugate Gradient method should be preferred over GMRES, since the coefficient matrix is a positive definite symmetric matrix. Tables 1 and 2, respectively, show iterations to tolerance (IT) and CPU times of algorithms V,  V , and  V5 described in Section 5 compared with: •  VL(geo) , which is the 2D Laplacian preconditioner introduced in Reference 14 inverted through one iteration of V-cycle with the geometric approach and Jacobi weight = 0.75 (as in Reference 14).   •  VL(gal) , which is the same as preconditioner  VL(geo) , but implemented through the Galerkin approach.

TA B L E 2 CPU times of the V-cycles V, V(̃), and the preconditioned GMRES with preconditioners
• V(̃) and  V(̃) , which are the same as V and  V but with Jacobi weight fixed as̃= 0.75 + √ min( , ) 4 (see Reference 14).
We do not consider any circulant preconditioner for two different reasons: first, in Reference 14 it has been shown that circulant preconditioners are slower than multigrid methods; second, it is well known that if used as preconditioners for multilevel Toeplitz matrices, multilevel circulant matrices cannot ensure a superlinear convergence character (see Reference 26).
In Tables 1 and 2, the numbers in bold highlight, in each row, the combination with the lowest computational time. We note that, as expected, when = , the convergence of V and  V is almost independent of the grid size and the amount of iterations is low. When , ≈ 0, the block-banded-banded-block preconditioner  V5 yields almost the same iterations as the full matrix  V , but with lower CPU times due to the lower computational cost per iteration. Moreover, the robustness of preconditioners  VL(geo) and  VL(gal) quickly deteriorates as , increase. Comparing V(̃) with V we note that the adaptive choice of the Jacobi weight explained in Section 5 allows slightly faster convergence with respect to the fixed weight̃.
When , ≈ 1, instead, the block-banded-banded-block preconditioner seems not to be suitable anymore. This is due to the decay of the coefficients of the discretized fractional operators which, as discussed in Section 5.2, tend to zero with order 2 − and 3 − , ∈ { , }, in the FVE and FV approaches, respectively. Notice that the slight improvement in both iterations and timings provided by the FV approach when , ≈ 1 is a direct consequence of the corresponding higher decay order. Tests not reported in Tables 1 and 2 show  V11 to be a robust solver, but still slower than V.
When ≠ , the number of iterations of all methods tends to increase as N x increases. This is due to the anisotropy of the diffusion along the two coordinate axes. Since hypothesis r s → c in Corollary 2 is not satisfied, neither is the approximation property, therefore the projectors in V-cycle should be built differently and a strategy like that proposed in Reference 25 should be explored. Nevertheless, using the GMRES with  V , not only halves the iteration with respect to V, but also seems to be much more robust in the anisotropic cases. Consequently, using the lighter preconditioner  V5 instead of  V allows to reach the lowest CPU times without losing in robustness. Now, let us fix N x = 2 11 − 1 and consider the solvers with the lowest CPU time in Table 2 for FV and FVE and for each combination of ( , ). More precisely, we consider solver  V5 for FV except for = = 0.5 and = = 0.9, where we use V, and solver V for FVE except for ( , ) = (0.3, 0.2), where we use  V5 . Figure 3 shows the 2-norm error versus the CPU time of such solvers for FV (solid line) and FVE (dashed line). We note that when , ≈ 0, the FV method is more efficient since it allows to compute solutions with smaller error than FVE in the same amount of time, despite the fact that for a given grid FVE is sometimes faster (see Table 2). FV seems to be more efficient than FVE in the anisotropic cases too, even for large , where FVE has a higher accuracy. Instead, in the isotropic cases with , ≈ 1 both approaches allow similar CPU times and, therefore, FVE becomes more suitable than FV.
We stress that due to the presence of the tridiagonal mass matrices, each matrix-vector product is more expensive in FVE than in FV. This goes in favor of FV since allows V-cycle to yield faster results than in the case of FVE, even when a larger number of iterations is required. Remark 9. Note that it is not possible to compare the iterations of  VL(geo) and V(̃) in Table 1 with preconditioners P 2,N and MGM 2D (J) in Reference 14, because therein the 2D discretization is different from the one given in Equation (16). Indeed, in References 13 and 14, the authors replaced the tridiagonal mass matrix with an identity matrix resulting in a mixed FV and FVE approach.

CONCLUSIONS AND FUTURE PERSPECTIVES
We have introduced a second-order FV method for problem (1) and we have numerically shown that it is a good alternative to the FVE approach when , ≈ 0. Moreover, we have proposed a block-banded-banded-block preconditioner for GMRES that allows a fast solution of the resulting linear systems in an amount of iterations to tolerance that is stable as the size of the coefficient matrix increases. When , ≈ 1, the FVE approach revealed more accurate than FV. In this case, a multigrid method used as standalone solver for the discretized problem should be preferred. Same as in Reference 14, we used damped Jacobi as smoother, but here we selected its weight adaptively, which yields better results if compared to the fixed weight proposed in Reference 14.
As highlighted in Remark 7, further improvements in terms of approximation error of the FV approach compared to the FVE approach when , ≈ 0 could be obtained using the shift p =  . The analysis of the related symbol and the design of symbol-based preconditioning and multigrid strategies will be subject of future studies.
which proves Equation (9) and concludes the proof. ▪

APPENDIX B. PROOF OF THEOREM 2
Proof of Theorem 2. Let us first show that F (x) is nonnegative, rewriting F (x) = t 1 (x)t 2 (x), with ) .