Conditioning of Finite Element Equations with Arbitrary Anisotropic Meshes

Bounds are developed for the condition number of the linear finite element equations of an anisotropic diffusion problem with arbitrary meshes. They depend on three factors. The first, factor proportional to a power of the number of mesh elements, represents the condition number of the linear finite element equations for the Laplacian operator on a uniform mesh. The other two factors arise from the mesh nonuniformity viewed in the Euclidean metric and in the metric defined by the diffusion matrix. The new bounds reveal that the conditioning of the finite element equations with adaptive anisotropic meshes is much better than what is commonly feared. Diagonal scaling for the linear system and its effects on the conditioning are also studied. It is shown that the Jacobi preconditioning, which is an optimal diagonal scaling for a symmetric positive definite sparse matrix, can eliminate the effects of mesh nonuniformity viewed in the Euclidean metric and reduce those effects of the mesh viewed in the metric defined by the diffusion matrix. Tight bounds on the extreme eigenvalues of the stiffness and mass matrices are obtained. Numerical examples are given.


Introduction
It has been amply demonstrated that significant improvements in accuracy can be gained when an appropriately chosen anisotropic mesh is used for the numerical solution of problems exhibiting anisotropic features.However, there exists a general concern in the scientific computing community that an anisotropic mesh, which can contain elements of large aspect ratio and small volume, may lead to ill-conditioned linear systems and this could outweigh the accuracy and efficiency improvements gained by anisotropic mesh adaptation.For isotropic mesh adaptation, Bank and Scott [2] (also see Brenner and Scott [3]) show that after proper diagonal scaling, the condition number of finite element equations with an adaptive mesh is essentially the same as that for a uniform mesh.Unfortunately, this result does not apply to anisotropic meshes nor to problems with anisotropic diffusion.
For problems with anisotropic diffusion and arbitrary meshes, several estimates have been developed for the extreme eigenvalues of the stiffness matrix.For example, Fried [7] shows that the largest eigenvalue of the stiffness matrix is bounded by the largest eigenvalues of element stiffness matrices.Shewchuk [18] obtains sharp bounds on largest eigenvalues of element stiffness matrices for linear triangular and tetrahedral finite elements.More recently, Du et al. [5] develop a bound that can be viewed as a generalization of Shewchuk's result to general dimensions and simplicial finite elements.
Estimation of the smallest eigenvalue for the general case appears to be more challenging.Standard estimates (e.g., see Ern and Guermond [6]) are linearly proportional to the volume of the smallest mesh element, which is typically too pessimistic for nonuniform meshes.Moreover, Apel [1,Sect. 4.3.3]shows that the order of the smallest eigenvalue of the stiffness matrix for a specific, specially designed anisotropic mesh is the same as for a uniform mesh.As a matter of fact, coefficient adaptive anisotropic meshes can even improve the conditioning for partial differential equations (PDEs) with anisotropic diffusion coefficients, as observed by D'Azavedo et al. [4] and Shewchuk [18,Sect. 3.2].A noticeable approach for obtaining sharper bounds for the smallest eigenvalue is proposed by Fried [7].The approach employs a continuous generalized eigenvalue problem with an auxiliary density function and its key is to find a lower bound for the smallest eigenvalue of the continuous problem.Bounds for the smallest eigenvalue of the stiffness matrix obtained with Fried's approach are valid for general meshes in any dimension but in d ≥ 3 dimensions they are less sharp than those obtained in this paper.
The objective of this paper is threefold.First, we develop tight bounds on the extreme eigenvalues and the condition number of the stiffness matrix for a general diffusion problem with an arbitrary anisotropic mesh.No assumption on the shape or size of mesh elements is made in the development.Our upper bound on the largest eigenvalue can be also expressed in terms of mesh nonuniformity viewed in the metric tensor defined by the diffusion matrix (which will hereafter be referred to as the mesh D-nonuniformity).It is comparable to those of Shewchuk [18] and Du et al. [5] but is expressed as a sum of patch-wise terms instead of element-wise terms as in the aforementioned references.The patch-wise nature gives a sharper bound and makes it more convenient to use in the development of diagonal scaling preconditioners.To obtain lower bounds for the smallest eigenvalue of the stiffness matrix we extend Bank and Scott's result [2] to arbitrary meshes.This generalization is not trivial and special effort has to be made to deal with the arbitrariness of the mesh.Along the way we establish anisotropic upper and lower bounds on the extreme eigenvalues of the mass matrix which are much tighter than estimates available in the literature.
The second objective of the paper is to provide a clear geometric interpretation for the obtained bounds on the condition number of the stiffness matrix.These bounds are shown to depend on three factors.The first factor is proportional to a power of the number of mesh elements and represents the condition number of the stiffness matrix for the linear finite element approximation of the Laplacian operator on a uniform mesh.The other two factors arise from the mesh nonuniformity in volume measured in the Euclidean metric (which will be referred to the mesh volume-nonuniformity) and from the mesh D-nonuniformity.
The third objective is to study diagonal scaling for the finite element linear system and its effects on the conditioning.We focus on the scaling with the diagonal entries of the matrix (Jacobi preconditioning) since it is an optimal diagonal scaling for a symmetric positive definite sparse matrix [10, Corollary 7.6 and the following].We show that the Jacobi preconditioning can eliminate the effects of the mesh volume-nonuniformity and improve those caused by the mesh D-nonuniformity, thus significantly reducing the effects of the mesh irregularity on the conditioning.From the practical point of view, this result indicates that a simple diagonal preconditioning can effectively transform the stiffness matrix into a matrix which has a comparable condition number as the one with a uniform mesh.
The outline of the paper is as follows.Section 2 briefly describes a linear finite element discretization of a general anisotropic diffusion problem.Estimation of the extreme eigenvalues and the condition number of the mass matrix is given in Sect.3. Section 4 deals with the estimation of the largest eigenvalue of the stiffness matrix.Bounds on the smallest eigenvalue and the condition number of the stiffness matrix and the effects of diagonal scaling are investigated in Sect. 5. A selection of examples in Sect.6 provides a numerical validation for the theoretical findings.Finally, conclusions and further remarks are given in Sect.7.

Linear finite element approximation
We consider the boundary value problem (BVP) of a general diffusion differential equation in the form where Ω is a simply connected polygonal or polyhedral domain in R d (d ≥ 1) and D = D(x) is the diffusion matrix.We assume that D is symmetric and positive definite and there exist two positive constants d min and d max such that where the less-than-or-equal sign means that the difference between the right-hand side and left-hand side terms is positive semidefinite.We are interested in the linear finite element solution of BVP (1).Assume that an affine family {T h } of simplicial decompositions of Ω is given and denote the associated linear finite element space by V h ⊂ H 1 0 (Ω).A linear finite element solution u h ∈ V h to BVP (1) is defined by Since both ∇u h and ∇v h are constant on K, we can rewrite the above equation as where D K is the integral average of D over K, i.e., In practice, the integrals in (3) and ( 4) have to be approximated numerically via a quadrature rule.Although this will change the definition of D K and the right-hand side term of (3) slightly, the procedure and the results in this paper will remain valid for this situation.Finite element equation ( 3) can be expressed in a matrix form.Denoting the numbers of elements and interior vertices of T h by N and N vi and assuming that the vertices are ordered in such a way that the first N vi vertices are the interior vertices, we have where φ j is the linear basis function associated with the j th vertex.In (5) and hereafter, we use the sum j with the index j ranging over all interior vertices, i.e., Substituting (5) into (3) and taking v h = φ i (i = 1, . . ., N vi ), we obtain the linear algebraic system where u = [u 1 , . . ., u N vi ] T and the stiffness matrix A and the right-hand side term f are given by and ∇φ i | K and ∇φ j | K denote the restriction of ∇φ i and ∇φ j on K. Our main goal is to estimate the condition number of stiffness matrix A.

Mass matrix
To start with, we consider the element mass matrix B for the reference element K (which is assumed to be unitary, i.e., | K| = 1), B = ( bi,j ), bi,j = K φi φj dξ, i, j = 1, . . ., d + 1, where φi 's are the linear basis functions associated with the vertices of K. Matrix B is symmetric and positive definite.Moreover, for the standard linear finite elements considered in this paper we have 1

Condition number of the mass matrix
Consider the global mass matrix Notice that where ω j is the element patch associated with the j th vertex and |ω j | is its volume.
The following theorem gives lower and upper bounds on the condition number of the mass matrix for any dimension and any mesh.Theorem 3.1 (Condition number of the mass matrix).The condition number of the mass matrix for the linear finite elements on a simplicial mesh is bounded by Proof.For an element K, let u K be the restriction of the vector u on K and B K the element mass matrix.Then, Rearranging the sum on the right-hand side according to the vertices and using (7), Similarly, we have

Relation to the estimates in the literature
If we denote the maximum number of mesh elements in a patch by p max , then and estimate (9) implies which is the bound obtained by Fried [7, inequality (24)].Moreover, for an isotropic mesh, where h max and h min are the largest and smallest element diameters.Substituting this into bound (10) gives which is precisely the standard estimate found in the literature (e.g., [6, Rem.9.10]).
For anisotropic meshes, on the other hand, the new estimate ( 9) is much tighter than both Fried's estimate (10) and the standard estimate (11) 9) also provides a tight lower bound, which is not available with (10) and (11).

Diagonal scaling for the mass matrix
It is known [10, Corollary 7.6 and the following] that for a symmetric positive definite sparse matrix, scaling by its diagonal entries (Jacobi preconditioning) is an optimal diagonal preconditioning (up to a constant depending on the maximum number of non-zeroes per column and row of the matrix).We are interested in a bound on the condition number after such preconditioning.
For a diagonal scaling S = (s j ), similarly to Theorem 3.1 we obtain and, for the Jacobi preconditioning s2 j = B jj , we have arrived at the following theorem by Wathen [19] who studies the effects of the diagonal scaling on the condition number of the Galerkin mass matrix.

Theorem 3.2 ([19, Table 1]). The condition number of the Jacobi preconditioned Galerkin mass matrix with a simplicial mesh has a mesh-independent bound
Theorems 3.1 and 3.2 show that the mesh volume-nonuniformity has a significant effect on the condition number of the mass matrix and this effect is completely eliminated by the Jacobi preconditioning.
As we will see later in Sect.5, diagonal scaling plays a similar role in reducing the effects of mesh nonuniformity on the condition number of the stiffness matrix.

Largest eigenvalue of the stiffness matrix
The following lemma is valid for any dimension.

Lemma 4.1 (Largest eigenvalue). The largest eigenvalue of the stiffness matrix A = (A ij ) for the linear finite element approximation of BVP (1) is bounded by
The largest eigenvalue of the diagonally (Jacobi) preconditioned stiffness matrix Proof.First, recall that for any symmetric positive semidefinite matrix M , Then, using the local indices on K and the definition of A jj from ( 6) and rearranging the sum according to the vertices, we have On the other hand, using the canonical basis vectors e j we have and altogether we get (12).Using the same procedure for a diagonal scaling S = (s j ) we obtain max For the Jacobi preconditioning we have s 2 j = A jj , which gives estimate (13).
Remark 4.2.Bound (13) can also be obtained by using the unassembled form of A as shown in [19,Sect. 3].However, the analysis employed in [19] cannot provide a lower bound on λ min (S −1 AS −1 ) other than the trivial one, λ min (S −1 AS −1 ) ≥ 0.

Geometric interpretation
Although Lemma 4.1 gives a very tight bound on λ max (A), it does not provide any explanation on how the mesh or the diffusion matrix affect the conditioning.We now derive a bound on λ max (A) in terms of mesh quantities and the diffusion matrix.Let F K : K → K be the affine mapping from the reference element K to the mesh element K, F K the Jacobian matrix of F K , j K the local index of φ j on K and φj K = F K • φ j K the corresponding basis function on K.
Using the chain rule, we have where Combining this with Lemma 4.1 yields Remark 4.3.If we denote the maximum number of elements meeting at a mesh point by p max , then bound (15) implies which is comparable to the estimates mostly found in the literature (e.g., [5,7,18]).Note that both this bound and (15) are less tight than the bound (12).
Remark 4.4.Bound (15) implies that the scaling will also lead to bounds similar to (13) and those in Sect. 5.In general, this scaling is greater than the Jacobi preconditioning (cf.( 14)) although they are equal in 1D or for a mesh that is uniform in the metric specified by D −1 (cf.Sect.4.1.3).
The ratio hK h min,K is a measure of the aspect ratio of K. Thus, for the case of D = I, the largest eigenvalue of A is bounded by the maximum volume-weighted element aspect ratio of the mesh.This is consistent with the observation by Shewchuk in [18] where a detailed discussion on the relation between the largest eigenvalue of the stiffness matrix and the element aspect ratio is available for the case of D = I in d = 2 and d = 3 dimensions.
For the general case D = I, on the other hand, it is more convenient to interpret bound (15) in terms of the mesh quality measures introduced in [11].We now proceed with this.

Mesh quality measures
The first measure is the alignment quality measure, which can be simply viewed as an equivalent to the aspect ratio of K in the metric specified by D −1 K .It is defined as and measures how closely the principal directions of the circumscribed ellipsoid of K are aligned with the eigenvectors of D K and the semi-lengths of the principal axes are proportional to the eigenvalues [15].Notice that In particular, Q ali,D −1 (K) = 1 implies that K is equilateral in the metric D −1 K .The second measure is the equidistribution quality measure defined as the ratio of the average element volume to the volume of K, both measured in the metric specified by where The equidistribution quality measure satisfies Notice that as the mesh is being refined.As a consequence, σ h can be considered as a constant.

Geometric interpretation (general case)
Using the quality measures we can rewrite the key factor and therefore Thus, λ max (A) is bounded by the maximum volume-weighted, combined alignment and equidistribution measure of the mesh in the metric D −1 K .When a mesh is adapted to the coefficients of the BVP, i.e., it is uniform in the metric D −1 , it will have the properties and Moreover, bound (18) will reduce to

Smallest eigenvalue and condition number of the stiffness matrix
The approach employed in this section was originally developed by Bank and Scott [2] for isotropic meshes.We generalize it here to arbitrary anisotropic meshes.
Hereafter, we will use C as a generic constant which can have different values at different appearances but is independent of the mesh, the number of mesh elements, and the solution of the BVP.
We start with bounds on λ min (A).
Lemma 5.1 (Smallest eigenvalue).The smallest eigenvalue of the stiffness matrix for the linear finite element approximation of BVP (1) is bounded from below by where | K| = 1 N |Ω| denotes the average element size.The smallest eigenvalue of the diagonally (Jacobi) preconditioned stiffness matrix is bounded from below by and Proof.Since Sobolev's inequality is different for d = 1, d = 2 and d ≥ 3 dimensions [8, Theorem 7.10], we treat these cases separately.
Case d = 1: Let C S be the constant associated with Sobolev's inequality.Using the inequality (2), Sobolev's inequality, and the equivalence of the vector norms, With scaling, In 1D, ∇φ j K | K = |K| −1 and therefore Using this in (24) gives (22).
Case d = 2: Consider a set of not-all-zero non-negative numbers {α K , K ∈ T h } (to be determined later) and a finite number q > 2. Let C P , C S , and C K be the constants associated with Poincaré's inequality, Sobolev's inequality, and the norm equivalence on K, respectively.Using (2), Poincaré's, Sobolev's and Hölder's inequalities and the norm equivalence for ûh , we have and therefore The largest lower bound on (25) is obtained for q = max 2, ln (N |K min |) with .
The choice q = 2 is viewed as the limiting case as q → 2 + .Estimate (21) follows from this, (25) and the definition of the average element size.
With scaling, we have For the Jacobi preconditioning s 2 j = A jj = K∈ω j ∇φ j • D K ∇φ j we choose where C φ = max i K =1,...,d+1 ∇ φi K 2 .With these and choosing the value for the index q in a similar manner as for the case without scaling we obtain (23).
Case d ≥ 3: This case is very similar to case d = 2. Again, from (2), Poincaré's, Sobolev's and Hölder's inequalities and the norm equivalence for ûh , we have Estimate (21) follows from this and the definition of the average element size.
The bound for the scaled stiffness matrix is obtained by choosing Combining Lemma 4.1, estimate (15) and Lemma 5.1 we obtain upper bounds on the condition number of the stiffness matrix and the scaled stiffness matrix.

Theorem 5.2 (Condition number of the stiffness matrix). The condition number of the stiffness matrix for the linear finite element approximation of BVP (1) is bounded by
and The condition number of the diagonally (Jacobi) preconditioned stiffness matrix is bounded by and 1, for d ≥ 3. (29)

Geometric interpretation
We now study the geometric interpretation of the bounds for the condition number.

Without scaling
Bounds ( 26) and ( 27) contain three factors, a base bound CN 2 d , a factor reflecting the effects of the mesh nonuniformity measured in the metric D −1 (mesh D-nonuniformity), and, if d ≥ 2, a factor reflecting the effects of the mesh nonuniformity in volume measured in the Euclidean metric (volume-nonuniformity).
The first factor N 2 d corresponds to the condition number of the stiffness matrix for the Laplacian operator on a uniform mesh (cf.Special Case 5.3 below).
The second factor reflects the effects of the mesh D-nonuniformity and can be understood as a volume-weighted, combined alignment and equidistribution quality measure of the mesh with respect to D −1 (cf.Sect.4.1.3).The third factor in ( 27) is It measures the effects of the mesh volume-nonuniformity (measured in the Euclidean metric) on the condition number.Notice that there is no effect in 1D and in 2D it is minimal.In d ≥ 3 dimensions the factor is proportional to the average of |K| −1+ 2 d over all elements.This is a significant improvement in comparison with previously available estimates which are proportional to

With scaling
Bounds (28) and (29) for the scaled stiffness matrix have the same base bound as without scaling.Hence, diagonal scaling has no effect on the condition number when the mesh is uniform and D = I.
Unlike (27), bounds ( 28) and ( 29) do not have the third factor which involves only the element volume (in comparison to the second factor which couples (F K ) −1 with D −1 ).In this sense, a properly chosen diagonal scaling can eliminate the effects of the mesh volume-nonuniformity on the condition number.Moreover, scaling can also significantly reduce the effects of the mesh Dnonuniformity.Indeed, the factors in ( 28) and (29) that couple (F K ) −1 with D −1 are asymptotically the L max{1, d 2 } (Ω) norm of (F K ) −1 D K (F K ) −T 2 whereas the corresponding factors in ( 26) and ( 27) are basically the maximum norm.
Furthermore, the D-related factor in (29) for d ≥ 2 can be rewritten in terms of the alignment quality measure Q ali,D −1 from Sect.4.1.2as Thus, the dependence of this D-related factor on the element volume is also mild: both Q ali,D −1 (K) and D K (the average of D over K) are invariant under the scaling transformation of K.
The following special cases are instructional to understand the interplay of the factors for different types of meshes.

Special Case 5.3 (Uniform meshes)
and bound (29) reduces to which is precisely the result of Bank and Scott [2, Theorems 4.2 and 5.2].In this case, the diagonal scaling becomes , where h j denotes the average length of the elements around the j th vertex.This scaling is equivalent to the change of basis functions Special Case 5.5 (Uniform meshes with respect to D −1 ).For a mesh that is uniform with respect to D −1 , i.e., coefficient adaptive, we have properties (19) and (20).Bounds (26)-(29) reduce to where σ h is defined in (17) and corresponds to the volume of the domain in the metric specified by D −1 .Thus, the condition number of the scaled stiffness matrix for a coefficient adaptive mesh has the optimal order of O(N Special Case 5.6 (Aligned meshes, d ≥ 2).For meshes aligned with the diffusion matrix but not necessarily fully coefficient adaptive (i.e., isotropic but not uniform with respect to D −1 ) we have From (30), bound (29) becomes Aside from the term depending on det(D), this bound is equivalent to (31).Hence, the diagonal scaling almost eliminates the effects of the mesh on the condition number for D-aligned meshes.
Special Case 5.7 (General M -uniform meshes).Finally, let us consider general M -uniform meshes, i.e., meshes that are uniform in the metric specified by a given metric tensor M which does not necessarily correspond to D −1 .In the context of mesh adaptation, an adaptive mesh is typically generated based on some estimate of the solution error and the associated metric tensor M is solution dependent.Thus, it is of interest to know what the impact of a given M on the conditioning of the stiffness matrix is.Recall [13] that an M -uniform mesh satisfies where M K is some average of M on K and σ h,M is defined as in (17) but with D replaced by M −1 .We have Hence, the bound on the condition number after diagonal scaling for an M -uniform mesh depends only on the volume-weighted average of For many problems such as those having boundary layers and shock waves, mesh elements are typically concentrated in a small portion of the physical domain.In that situation, we would expect that M differs significantly from D −1 only in small regions.As a consequence, the volume-weighted average of over the whole domain may remain small and therefore the condition number of the scaled stiffness matrix for anisotropic adaptive meshes does not necessarily increase as much as generally feared.
This effect can be observed in Examples 6.2 and 6.3.Figures 3a and 4a show that the effects of anisotropic adaptation are completely neutralized by the diagonal scaling when the number of anisotropic elements is small in comparison to N .

Numerical experiments
In this section we present numerical results for a selection of one-, two-, and three-dimensional examples to illustrate our theoretical findings.
Note that all bounds on the smallest eigenvalue contain a constant C. We obtain its value by calibrating the bound for λ min (S −1 AS −1 ) with Delaunay (Example 6.4) or uniform meshes (all other examples) through comparing the exact and estimated values.For the largest eigenvalue we use explicit bounds (12) and (13).
First, we give examples with predefined meshes to demonstrate the influence of the number and shape of mesh elements on the condition number of the stiffness matrix and to verify the improvement achieved with the diagonal scaling.For the tests, we employ the Laplace operator (i.e.D = I) and a mesh on the unit interval, square, and cube, for 1D, 2D, and 3D, respectively.Example 6.1 (d = 1, D = I, Chebyshev nodes).For a simple one-dimensional example we choose a mesh given by Chebyshev nodes in the interval [0, 1], The exact condition number of the stiffness matrix and its estimates (26) and (28) are shown in Figs.1a (without scaling) and 1b (with scaling) while those for the extreme eigenvalues and their estimates are given in Figs.1c (without scaling) and 1d (with scaling).Figure 1a shows that the estimate (26) is much sharper than the standard estimate with λ min (A) ∝ |K min |.The former has the same asymptotic order as the exact value as N increases, whereas the latter is too pessimistic and has a higher asymptotic order.The difference is caused by the estimate of the smallest eigenvalue (Fig. 1c).Notice that the estimates on the largest eigenvalue are very tight, both for the scaled and the unscaled cases.
The results clearly show the benefits of diagonal scaling: the order for the condition number of the scaled stiffness matrix in Fig. 1b is O(N 2 ln N ), which is almost the same as for uniform meshes, whereas that without scaling in Fig. 1a is O(N 3 ).It can be shown analytically that the orders of the nonuniformity factors in (26) and (28) for the Chebyshev nodes defined with (32) are O(N ) and O(ln N ) and those of the corresponding condition numbers are O(N 3 ) and O(N 2 ln N ).
Thus, the numerical and theoretical results are consistent and the improvement by diagonal scaling from the maximum norm to the L 2 norm is significant in this example.Example 6.2 (d = 2, D = I, anisotropic elements in a unit square).For this 2D example we use a mesh for the unit square [0, 1] × [0, 1] with O(N 1/2 ) skew elements, as shown in Fig. 2a.First, we fix the maximum aspect ratio at 125 : 1 and increase N to verify the dependence of the condition number on N (Fig. 3a).Then, we fix N at 20,000 and change the maximum aspect ratio of the mesh elements to investigate the dependence of the conditioning on the mesh shape (Fig. 3b).
Figure 3a shows the averaging effect of the diagonal scaling: the scaling significantly reduces the condition number and, when N becomes large enough, the conditioning of a scaled system is     comparable to the condition number on a uniform mesh.Moreover, the estimated value of the condition number with or without scaling has the same order as the exact value as N increases.
Figure 3b provides a good numerical validation of (27), namely that the condition number of the unscaled stiffness matrix is linearly proportional to the largest aspect ratio2 .With scaling, the condition number is still increasing with an increasing aspect ratio, since the average aspect ratio (in accordance to (29)) is also increasing.Nevertheless, the condition number after scaling is smaller by a factor of 10.
Figure 3b also shows that our estimate of the condition number with scaling has the same (linear) order as the exact value as the maximum aspect ratio increases, whereas the bounds for the unscaled case has a slightly higher order.This indicates that the estimation can be further improved.
As for the estimates on the extreme eigenvalues, the results are mainly the same as in Example 6.1.For this reason, we omit them in 2D and 3D to save space.Example 6.3 (d = 3; anisotropic elements in a unit cube).In this example, we repeat the same test setting as in Example 6.2: fixed anisotropy (25 : 25 : 1) with increasing number of elements (Fig. 4a) and a fixed N = 29, 478 paired with the changing anisotropy of the mesh (Fig. 4b).The results shown in Fig. 4 are essentially the same as in 2D.Since the mesh used in this example has        3: Condition number before and after scaling for a predefined 3D mesh (Fig. 2b) as a function of (a) the number of mesh elements and (b) the maximum element aspect ratio.
better than in the case of a quasi-uniform mesh (Fig. 6a), confirming the observation that, depending on the problem, a quasi-uniform mesh is not necessarily the best mesh from the conditioning point of view.For both cases, diagonal scaling does not improve the condition number significantly.This is expected since both meshes are almost volume-uniform.To explain why the mesh in Fig. 6b is (almost) volume-uniform, we recall from ( 16) and ( 19) that an M -uniform mesh with respect to The diffusion matrix D in (33) satisfies det(D) = 1000.Thus, |K| = const.
The largest condition number is in the case of the purely solution-adaptive mesh (Fig. 6c).This is because the mesh is not volume-uniform and its elements are not aligned with D −1 .Since the mesh is far from being uniform in size, scaling will have a significant impact, as it can be verified in Fig. 6c: the condition number after the scaling is even smaller than the condition number with Delauney meshes.
Conditioning with a mesh that is both coefficient-and solution-adaptive (Fig. 6d) is not as good as in the case of the purely coefficient-adaptive mesh but better than in the case of the purely adaptive and Delaunay meshes.
In all four cases we observe that the developed estimates for the condition number of the stiffness matrix are reasonably tight and have the same order as the exact values as N increases for both unscaled and scaled cases.

Mass matrix
Our new estimate (8) of the condition number of the Galerkin mass matrix is tight within a factor of (d + 2) from both above and below for any mesh with no assumptions on mesh regularity or topology.It this sense, it is optimal and truly anisotropic.

Stiffness matrix
Lemma 4.1 provides an estimate of the largest eigenvalue of the stiffness matrix which is simple to compute and is tight within a factor of (d + 1) from both above and below for any mesh.This is in contrast to many existing estimates which are proportional to the maximal number of elements meeting at a mesh point.
New bounds (21)-(23) on the smallest eigenvalue and ( 26)-(29) on the condition number of the stiffness matrix are a significant improvement in comparison to the previously available estimates.
First, the new bounds show that the conditioning of the stiffness matrix with an arbitrary (anisotropic) mesh is much better than generally assumed, especially for d = 1 and d = 2.
Second, the new bounds are truly anisotropic and valid for any mesh since no assumptions on the mesh regularity were made.
Third, bounds (26) and ( 27) reveal what affects the conditioning.There are three factors.The first (base) factor CN 2 d describes the direct dependence of the condition number on the number of mesh elements and corresponds to the condition number for the Laplace operator on a uniform mesh.The second factor describes the effects of the mesh D −1 -nonuniformity, i.e., the interplay between the shape and size of mesh elements and the coefficients of the BVP.It is O(1) for a coefficient-adaptive mesh, i.e., a mesh satisfying (19).The third factor measures how the mesh volume-nonuniformity further affects the condition number.It has no effect in 1D, a minimal one in 2D, and a substantial effect in 3D and higher dimensions.This means that even if the mesh is coefficient-adaptive and the second factor is O(1), the mesh volume-nonuniformity can still have a significant impact on the condition number for d ≥ 3.
Fourth, a simple diagonal scaling, such as the Jacobi preconditioning, can significantly improve the conditioning.Bound (29) for the condition number after scaling does not contain the factor for the mesh volume-nonuniformity.As a consequence, for a coefficient-adaptive mesh, this bound reduces to the base factor CN 2 d .In this sense, diagonal scaling eliminates the effects of the mesh volume-nonuniformity.It can also significantly reduce the effects of the mesh nonuniformity with respect to D −1 : the influence reduces essentially from the maximum norm to the L max{1, d 2 } norm of (F K ) −1 D K (F K ) −T 2 .Moreover, for a preconditioner that is invariant to diagonal scaling it follows that the condition number of the preconditioned stiffness matrix is typically smaller than κ(S −1 AS −1 ) which in turn has a much lower bound than κ(A) (cf.( 27) and (29)).For example, consider an incomplete Cholesky decomposition of A, A = LL T + E.

It follows that
is actually an incomplete Cholesky decomposition of S −1 AS −1 since S −1 L has the same sparsity pattern as L. Then from the identity we see that the preconditioned matrix of A with preconditioner L is equivalent to the preconditioned matrix of S −1 AS −1 with preconditioner S −1 L. As a result, the performance of the preconditioning technique on A is the same as that on S −1 AS −1 which has a much smaller condition number than A.
Although there is no estimate yet on the condition number of the preconditioned system, the above observation may provide a partial explanation for the good performance of ILU preconditioners with anisotropic meshes observed in [12].Numerical experiments (Figs.3a and 4a) indicate that although the new bounds have the same order as the exact value as the number of elements increases, they may have higher asymptotic orders than the exact value as the element aspect ratio increases.These may deserve further investigations.
Finally, we would like to point out that although the study in this paper has been done specifically for the linear finite element discretization, the approach can be generalized for higher order finite elements without major modifications.

4. 1 . 1
Special case D = I For the simplest case of D = I, bound (15) has a rather simple interpretation.The quantity (F K ) −1 2 can be bounded by the reciprocal of the in-diameter h min,K of K [15, Lemma 5.1.2].If we denote the average aspect of K by hK (i.e., hK = |K| 1 d ), then we can rewrite (15) as For a uniform mesh and D = I, bounds (26)-(29) yield κ(A) ≤ CN 2 d and κ(S −1 AS −1 ) ≤ CN 2 d , which is the base bound.Hence, the diagonal scaling has no effect on the condition number when the mesh is uniform and D = I.Special Case 5.4 (Isotropic meshes, D = I, d ≥ 2).For an isotropic mesh and D = I, ) estimate (standard) κ(A) estimate κ(A) exact (a) Condition number.

Figure 1 :
Figure 1: Example 6.1: Exact and estimated condition number and eigenvalues of the stiffness matrix as a function of N (d = 1).

Figure 3 :
Figure 3: Example 6.2: Condition number before and after scaling for a predefined 2D mesh (Fig. 2a) as a function of (a) the number of mesh elements and (b) the maximum element aspect ratio.

Figure 4 :
Figure 4: Example 6.3: Condition number before and after scaling for a predefined 3D mesh (Fig.2b) as a function of (a) the number of mesh elements and (b) the maximum element aspect ratio.