# Question

When it is about solving a system of equations, is it possible to predict that whether high-frequency noise (e.g. checker-boarding) is likely to appear in the converged solution by looking at coefficient matrix attributes (e.g. definiteness, etc)? What solver types are less likely to result oscillations, direct solvers or iterative solvers?

# Description

There is a stress analysis problem where the momentum equation is discretised by a coupled finite volume method. The coefficient matrix is not positive definite. Oscillations (checker-boarding) appear in the final solution when tetrahedral meshes are used. Bi-conjugate gradient and GMRes are used to solve the system. GMRes is much slower than the former.

• Comments are not for extended discussion; this conversation has been moved to chat. – nicoguaro Aug 6 '20 at 15:23

Most of this was already discussed in the comments, but I would like to elaborate and put a detailed answer.

There are no elementary characteristics (definiteness, symmetry, bandwidth) which can tell you whether the underlying (mixed or not) FEM/FVM is stable to solve the continuous problem. You can not tell anything about that just by looking at those characteristics. You should look at consistency, coercitivity (and if applies, inf-sup stability too) and boundedness of the corresponding finite dimensional operator.

There are some simple cases though. For example, if the matrix SPD then it is non-singular and in general, you can replace coercivity of the FEM with positive definiteness of the matrix. However, testing for positive definiteness is no easier than solving the linear system -Cholesky factorization, SVD, or inverse power iteration to find the minimum eigenvalue all are O(n^3)-.

However, you can also find many counter-examples which have zero eigenvalues but do not result in any spurious oscillations. For example,

1. Take the Poisson problem in its primitive form with all Neumann boundary conditions. In this case, you will find a solution but it will unique up to an additive constant, i.e. if $$u$$ is a solution then $$u+c$$ is a solution for any $$c\in\mathbb{R}$$. This because there is one zero eigenvalue and the null space contains the constant element. Another similar case is enclosed flow problems, where the pressure is unique up to an additive constant and the coefficient matrix has a zero eigenvalue. So the existence of a zero eigenvalue is not enough of say there will be spurious oscillations, but if there are no zero eigenvalues then it is a good indicator.

2. You can not check whether the eigenvalues of the coefficient matrix are going to zero as you refine the mesh as well, because that is what happens for many problems. Take $$P_2-P_1$$ discretization of the Stokes equations -note that, it is inf-sup stable-, if you compute the eigenvalues as you refine the mesh you will see that the minimum eigenvalue will tend to zero. In this case, you will not see any spurious oscillations. So the convergence to zero eigenvalue is not enough to say that there will be spurious oscillations.

Now, there were some objections to this answer in the comments, which are fair questions. First one was the inf-sup condition for the mixed methods. Numerically testing it is possible as pointed by Chenna K ("The inf-sup" test by Bathe and Chapelle, Computers & Structures, 47(4–5):537-545, 1993.). Let's understand what this paper says but before that:

You should know whether the method you are using is inf-sup stable before you use it! You should not use theoretically unproven techniques to solve problems. If you are doing research and you want to spend computer time rather than your own to eliminate bad approaches, it is fine. However, for the sake of everything nice and good, do not use theoretically unproven methods to solve real life problems. At best, you will make bad predictions and waste time spent doing experiments, at worst, you will kill people because of your bad predictions. Take many, many failed rocket launches because of bad numerics as a lesson.

Let's get back to the paper. The inf-sup condition -or as sometimes called Ladyzhenskaya-Babushka-Brezzi (LBB) condition- is a sufficient condition for saddle point problems to have a unique solution. It is criteria for stability of saddle point problems. Mixed methods give rise to saddle point problems. So let's discuss what it is. Say we are interested in solving the following problem: Find $$u\in V$$ and $$p\in Q$$ (for suitable function spaces $$V$$ and $$Q$$), such that \begin{align} &a(u,v) + b(p,v) &= L_1(f) \\ &b(u,q) &=L_2(g) \end{align} for all $$v\in V$$ and $$q\in Q$$. We need $$a(\cdot,\cdot)$$ to be coercive (there are some details to this, maybe check Boffi, Brezzi, Fortin "Mixed Finite Element Methods and Applications"), along with continuity of $$a(\cdot,\cdot)$$ and $$b(\cdot,\cdot)$$, then if $$b(\cdot,\cdot)$$ satisfies the inf-sup condition \begin{align} \inf_{q\in Q}\sup_{v\in V, v\neq 0} \tfrac{b(v,q)}{||v||_V||q||_Q} \geq \beta > 0 \end{align} then there is a unique solution pair to the saddle point problem. $$||\cdot||_X$$ represent the endowed norm of the function space $$X$$. This is, of course, continuous and it is impossible to compute this numerically. However, there is an equivalent discrete problem we can solve. Let $$W$$ be a matrix such that $$||u_h||_V= (U_h, WU_h)$$ for all $$u_h\in V_h$$ where $$V_h$$ is a finite dimensional subspace of $$V$$, $$u_h = \sum_i U_h(i)\phi_i$$ with $$U_h(i)\in\mathbb{R}$$ and $$\phi_i$$ are the elements of the basis of $$V_h$$. Define $$R$$ a matrix similarly, $$||p_h||_V= (P_h, RP_h)$$ for all $$p_h\in Q_h$$. Now you can write the inf-sup condition above in discrete form \begin{align} \inf_{q_h\in Q_h}\sup_{v_h\in V_h, v_h\neq 0} \tfrac{b(v,q)}{||v||_V||q||_Q} \geq \beta > 0 \end{align} or \begin{align} \inf_{P_h}\sup_{U_h, ||U_h||\neq 0} \tfrac{(BU_h,P_h)}{(U_h, WU_h)(P_h, RP_h)} \geq \beta > 0. \end{align} $$\beta$$ is a positive real constant independent of the mesh size, but may depend on other qualities of the mesh, e.g. aspect ratio; $$B$$ is a matrix defined similarly to $$W$$ and $$R$$ but represents $$b(\cdot,\cdot)$$. Still not that solvable, can be restated as an eigenvalue problem: Find minimum eigenvalue $$\lambda_{min}$$ different than zero which satisfies $$BW^{-1}B^TQ_h = \lambda RQ_h,$$ then $$\beta=\lambda_{min}$$. This is proposition 2.2 in the aforementioned paper, see proposition 2.1 for an equivalent generalized eigenvalue problem. Now, we have a numerical test to see if a given FE pair is inf-sup stable. I will argue that although a good tool, it is not really an answer to your question.

1. You need to solve a generalized eigenvalue problem that is barely related to the original coefficient matrix. If $$V$$ is endowed with an energy norm, maybe $$W=A$$ where $$A$$ is the matrix corresponding to the bilinear form $$a(\cdot,\cdot)$$ and $$R$$ is the mass matrix defined on $$Q_h$$, but this information is not available from the coefficient matrix.

2. If you can show that $$\lambda_{min}$$ does not change as you uniformly refine the mesh, then you know that the discretization is inf-sup stable. You are all done right? No, solving this eigenvalue problem still requires a lot of computational work -equivalent of solving the original linear system- and you don't have immediate access to these matrices which can be full (see the discussion right before subsection 2.3). The matrix size will grow as you refine the mesh, so this test is very limited in its application.

3. Say you have the resources to do the computation and you found that $$\lambda_{min}$$ seems like it is converging to zero. How do you know if it is not due to round-off errors or underflow or catastrophic cancellation? You can increase the precision and check again I guess, but it is becoming really ridiculous now.

Summary of this section: Just check the theoretical results, if the method is stable theoretically you are good. You can do numerical tests, which are pretty cool honestly, but they are also quite a bit work and only the affirming results are reliable.

Now, Chenna K also brought up the issue of inappropriate choice of numerical integration. This is a variational crime, it will introduce some extra error in the solution. It might happen that the quadrature error dominates and corrupt the solution resulting in spurious oscillations. However, you can't tell that just by looking at the coefficient matrix. Multiple reasons

1. Given a "correct" reference matrix, I can give you two other matrices of the same size; one with an inappropriate numerical integration scheme (either changing the order or changing the quadrature points and weights) and one with different basis functions $$\psi_i$$ for the same finite element space, and you would not be able to tell which is which just by comparing these matrices to the reference matrix. At least, not for cheap.

2. The coefficient matrix does not have to change. You can exactly compute the coefficient matrix, but if you don't compute the right hand side vector sufficiently accurately then the solution may still be wrong and/or contain spurious oscillations.

In summary of this part; you can not figure out if the integration scheme is inappropriate just by looking at the coefficient matrix.

In summary of this whole answer: There is no property of a matrix which will inform you about stability of the underlying discretization. There are some numerical tests but they are no cheaper than solving the problem and you can not use them to reject a method without putting significant computational effort. You can maybe use those tests to affirm the stability of a method, but that affirmation should be later proven by theory.

• "Existence of a zero eigenvalue is not enough of say there will be spurious oscillations", but if there are no zero eigenvalues is it true to say that the solution will be unique? – Alish Aug 6 '20 at 9:16
• It will be unique but not necessarily oscillation free. Take pressure robustness, for example. For some mixed finite element methods used for the solution of Navier-Stokes equations, the error in the solution of velocity field depends on inverse of viscosity times the error in the pressure. If the viscosity is small enough or the pressure error is large enough, there will be no zero eigenvalues, the solution to the linear system will be unique but the velocity solution may be garbage. – Abdullah Ali Sivas Aug 6 '20 at 15:20