Study on Matrix Analysis For Steady Problems

On a triangular mesh under the assumption that r The results obtained in [63, 64, 339] have illustrated the significance of DMPs for the analysis and design of numerical approximations. Various extensions and generalizations were published during the past three decades, see [61, 100, 151, 179, 180, 341, 311] and references therein. The frequently cited monograph by Ikeda [167] is devoted entirely to DMP for finite element models of convection-diffusion phenomena. Some low-order approximations of transport equations are known to satisfy a DMP unconditionally or under rather mild restrictions on the angles or aspect ratios of mesh cells. However, most a priori proofs are based on a set of sufficient conditions which become overly restrictive in the case of higher-order discretizations, singularly perturbed convection-diffusion equations, and anisotropic diffusion problems. A possible remedy to this problem is proposed in the next chapter. In this section, we review the algebraic constraints that ensure DMP and/or positivity preservation for steady transport problems of elliptic and hyperbolic type. A brief summary of the corresponding geometric conditions will also be presented.

A key ingredient of the mathematical theory behind the discrete maximum principles and positivity preservation is the following monotonicity concept [63, 339]. Definition 3.8. A regular matrix A = {ai j} is called monotone if A−1 _ 0. This kind of monotonicity is equivalent to the requirement that, for any vector u, Au _ 0 ) u _ 0. Matrix Analysis for Steady Problems 109 Of course, it is impractical to compute the inverse of A and check the sign of its entries. Instead, we will deal with a special class of matrices which are known to be monotone under certain constraints on the sign and magnitude of their coefficients. Definition 3.9. A regular matrix A = {ai j} is called an M-matrix if A−1 _ 0 and ai j _ 0, 8 j 6= i.

Available online at www.ignited.in Page 2

In other words, an M-matrix is a monotone matrix with nonpositive off-diagonal entries. These properties ensure positivity and convergence of iterative solvers. Definition 3.10. A matrix A = {ai j} is called diagonally dominant (by rows) if |aii j 6=i |ai j|, 8i. (3.53) Such a matrix is called strictly diagonally dominant if all inequalities are strict |aii j 6=i |ai j|, 8i. (3.54) Definition 3.11. A matrix A = {ai j} of size N×N is called irreducible if there is no N ×N permutation matrix P such that the following transformation is possible PAPT = _ A11 A12 0 A22 _ , where the size of A11 isM×M, the size of A22 is (N−M)×(N−M), and 1_M

Consider the steady transport-reaction equation (3.9) discretized by a finite difference, finite volume, or finite element scheme. Let the approximate solution uh, where the subscript h refers to the mesh size, be determined by a finite number ¯N of degrees of freedom u1, . . . ,u¯N that represent pointwise nodal values, control volume averages, or coefficients of piecewise-polynomial basis functions, respectively. Hence, all information about the solution uh can be packed into a vector u 2 R ¯N .

Lh : R

!R

Regardless of the underlying approximation technique, we define this mapping as Lhu = Au, where A = {ai j} is a sparse ¯N × ¯N matrix and u = {ui} is the vector of nodal values. The sparsity pattern of A depends on the mesh, on the type of discretization, and on the numbering of nodes. Since some nodal values are known from the Dirichlet boundary conditions, the size of the algebraic system reduces accordingly. Let the first N nodes be associated with the unknown degrees of freedom, and the rest with the Dirichlet boundary values. This numbering convention implies that the 108 3 Maximum Principles row/column numbers NN +1, . . . , ¯N }, respectively. Thus, uu1, . . . ,uN} is the vector of unknowns, whereas uuN+1, . . . ,u¯N } is given by the prescribed boundary data ug. (3.49) In this notation, the system of algebraic equations for the components of uAubAg, (3.50) where bNeumann boundary conditions, if any. In a practical implementation, it is convenient to incorporate the Dirichlet boundary conditions into the ¯N × ¯N matrix A and solve the extended linear system [63] In other words, the ¯N × ¯N matrix ¯A is constructed from A by setting AAI, where I denotes the identity matrix with ¯N −N rows and columns. If the solution of the continuous problem satisfies Theorem 3.5 or 3.7, it is natural to require that the maxima of ubounded by those of ug. Likewise, all nodal values should be nonnegative if Theorem 3.6 or 3.8 is applicable. To verify the validity of DMP, one needs to analyze the properties of the discrete operator ¯A. Remark 3.8. In the context of linear systems, irreducibility ensures that it is impossible to extract a subsystem that can be solved independently. Matrices that result from discretization of partial differential equations are irreducible in most cases. Definition 3.12. A matrix A = {ai j} is irreducibly diagonally dominant if it is irreducible and diagonally dominant, with strict dominance for at least one row i. The following theorem yields a set of sufficient conditions which are commonly employed in DMP analysis based on the M-matrix property of discrete operators.

110 3 Maximum Principles

Available online at www.ignited.in Page 3

Proof. A common approach to the proof of this theorem is based on the splitting A = D−C, where the diagonal part D = diag(A) > 0 is nonsingular and C _ 0. The diagonal dominance makes it possible to prove that the spectral B = D−1C B) < 1. This condition holds if and only if the series (I−B)−1 = I+B+B2+B3+. . .converges, see [131, 339] for technical details. Hence, A−1 = (I−B)−1D−1 _ 0. _ If all diagonal entries of A are strictly positive and there are no positive offdiagonal ones, then diagonal dominance (3.53) requires that all row sums be nonnegative. The following definition summarizes the corresponding sign conditions. Definition 3.13. A matrix A = {ai j} is said to be of nonnegative type [64, 69] if aii > 0, 8i, (3.55) ai j _ 0, 8 j 6= i, (3.56)

ai j _ 0, 8i. (3.57)

Note that conditions (3.55)–(3.56) impose the same constraints as the third basic rule from Section 1.6.3. The second basic rule is satisfied if (3.57) holds as equality. Under the assumptions of Corollary 3.12, the nonnegativity conditions are sufficient (but not necessary) for the matrix A to be monotone. Some other useful criteria related to M-matrices and monotonicity can be found in [46, 132, 153, 311, 346].

Given a discretization of the form (3.51), the monotonicity of the matrix ¯A makes it possible to prove discrete counterpart s of all maximum, minimum, and comparison principles established in Section 3.1. The uniqueness of the solution vector u follows from the regularity of ¯A. The usual approach to the proof of monotonicity is based on Theorem 3.13 and Corollary 3.12 since the nonnegativity conditions (3.55)–(3.57) are easy to verify for an arbitrary space discretization of the transport equation. To prove that the solution of problem (3.51) nodes, we need a discrete counterpart of the v _ 0. At the continuous level, it implies that Lu = L(u+c) for an arbitrary constant c. According to the second basic rule from Section 1.6.3, the discrete operator A should have zero row sums to inherit this property. Thus, the global discrete maximum principle for nodal values can be formulated as follows. 3.2 Matrix Analysis for Steady Problems 111

That is, if the source term is absent and the solution has a ui must also assume this value. This property is guaranteed by the zero row sum condition (3.58). Only a poor discretization of the transport equatv _ 0 would produce ui 6= ¯ u in this situation ([268], p. 39). Since estimate (3.63) holds for any interior node, successive application of the local DMP can be used to prove (3.59) if Asums of ¯A are nonvanishing for some i DMP may cease to hold but positivity preservation can be inferred from the fact that ¯A is monotone. If the matrix ¯A is given by (3.52), where Ais monotone and A0, then discretization (3.51) is positivity-preserving, that is, b _ 0 ) u _ 0. (3.64) Proof. The matrix ¯A given by (3.52) is regular if and only if the block AFurthermore, the inverse matrices. Some finite element approximations to elliptic problems like (3.47) are known to satisfy the DMP conditions (3.55)–(3.57) on a suitably designed mesh. The derivation of geometric constraints that ensure monotonicity has been one of the primary research directions in the DMP analysis for finite element schemes [64, 100, 179, 180]. Below, we present some useful geometric criteria in the context of linear and bilinear Galerkin discretizations of the Laplace operator in two space dimensions. Definition 3.14. A triangular mesh is called strongly acute weakly acute (or nonobtuse Theorem 3.21. The discrete Laplace operator ¯A for the linear finite element approximation on a triangular mesh of weakly acute type is monotone [18, 64]. This classical result dates back to the paper by Ciarlet and Raviart [64]. In the 3D case, a tetrahedral mesh is said to be of acute type if all internal angles between the faces of tetrahedra the discrete Laplace operator is monotone if linear finite elements are employed [189]. Definition 3.15. A triangular mesh is a Delaunay triangulation if no vertex of this mesh is inside the circumcircle of any triangle to which it does not belong.

Available online at www.ignited.in Page 4

Delaunay triangulations maximize the minimum angle as to avoid excessively stretched triangles. It is known that there exists a unique Delaunay triangulation for 114 3 Maximum Principles any set of points that do not lie on the same line.Moreover, fast algorithms are available for creating such triangulations [55, 115, 226], which makes them very popular with finite element practitioners. Figure 3.4 displays a simple 2D Delaunay triangulation generated using the MATLAB function delaunay. In three dimensions, no vertex of a tetrahedron is inside the circumsphere of any other tetrahedron. However, the discrete Laplacian operator for a linear FEM approximation on an arbitrary 3D Delaunay triangulation may fail to be a matrix of nonnegative type [18]. This does not necessarily cause a violation of the DMP but it cannot be ruled out anymore. Definition 3.16. A rectangular mesh is called nonnarrow if the ratio of longest and shortest mesh edge is not greater than p2 for any rectangle [100].

This theorem explains why iterative solution techniques that rely on theM-matrix property may experience convergence problems when applied to discretizations of second-order PDEs on quadrilateral/hexahedral meshes with high aspect ratios. Geometric DMP conditions for various elliptic and parabolic problems have been formulated building on the above results [100, 116, 179, 180]. Even if convective effects are present, it is desirable to use a sufficiently regular mesh that satisfies the above conditions, so that at least the discrete diffusion operator poses no hazard to DMP. Moreover, it may offset a nonmonotone convective part if the Peclet number is small or a sufficiently large amount of artificial diffusion is added [51, 69]. Alternatively, the upwind triangle method or other techniques can be used to construct a monotone approximation of convective terms [11, 167, 200, 288, 348].

1. M. Ainsworth and J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis. John Wiley & Sons, New York, 2000. 2. M. Ainsworth, J.Z. Zhu, A.W. Craig, O.C. Zienkiewicz, Analysis of the Zienkiewicz-Zhu aposteriori error estimator in the finite element method. Int. J. Numer. Methods Engrg. 28:9 (1989) 2161–2174. 3. J.D. Anderson, Jr., Computational Fluid Dynamics. The Basics with Applications. McGraw- Hill, 1995. 4. J.D. Anderson, Jr., Modern Compressible Flow: With Historical Perspective, McGraw-Hill, 1990. 5. F. Angrand, A. Derieux, L. Loth, G. Vijayasundaram, Simulation of Euler transonic flows by means of explicit finite element type schemes. INRIA Research Report 250, 1983. 6. F. Angrand and A. Derieux, Some explicit triangular finite element schemes for the Euler equations. Int. J. Numer. Methods Fluids 4 (1984) 749–764. 7. P. Arminjon and A. Dervieux, Construction of TVD-like artificial viscosities on 2-dimensional arbitrary FEM grids. INRIA Research Report 1111, 1989. 8. M. Arora and P.L. Roe, A well-behaved TVD limiter for high-resolution calculations of unsteady flow. J. Comput. Phys. 132 (1997) 3–11. 9. Athena test suite, http://www.astro.virginia.edu/VITA/ATHENA/dmr.html. 10. K. Baba and M. Tabata, On a conservative upwind finite element scheme for convective diffusion equations. RAIRO Numerical Analysis 15 (1981) 3–25. 11. I. Babuˇska and W.C. Rheinboldt, Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. 15:4 (1978) 736–354. 12. I. Babuˇska and W.C. Rheinboldt, A posteriori error estimates for the finite element method. Int. J. Numer. Methods Engrg. 12:10 (1978) 1597–1615. 13. W. Bangerth and R. Rannacher, Adaptive finite element methods for differential equations. Lectures in Mathematics, ETH Z¨urich, Birkh¨auser, 2003. 14. R.E. Bank, A.H. Sherman, A. Weiser, Some refinement algorithms and data structures for regular local mesh refinement. In: R. Stepleman (ed.), Scientific Computing, Applications of Mathematics and Computing to the Physical Sciences, IMACS Transactions on Scientific

Available online at www.ignited.in Page 5

Computation, Vol. I. North-Holland, Amsterdam, 1983, 3–17. 15. R.E. Bank and R.K. Smith, Mesh smoothing using a posteriori error estimates. SIAM J. Numer. Anal. 34 (1997) 979–997. 16. T.J. Barth, Numerical aspects of computing viscous high Reynolds number flows on unstructured meshes. AIAA Paper, 91-0721, 1991. 17. T.J. Barth, Aspects of unstructured grids and finite volume solvers for the Euler and Navier- Stokes equations. In: Lecture Series 1994-05, von Karman Institute for Fluid Dynamics, Brussels, 1994.