# Roundoff errors in FEM computations - generalized eigenvalues

This is a continuation of my previous question. I am trying to effectively compute a bound for the roundoff errors in some FEM computation (2d polygons, triangular meshes). Below I will write some of my computations which lead to a quantifiable bound (if I'm not mistaken), but the error is quite large for small mesh size $$h$$. Therefore, my questions are:

• is there a better way to quantify this type of errors? Do you know references who do this type of computation?

• is the error bound obtained below comparable to what is observed in practice? What is the general admitted threshold $$h_0$$ such that for $$h the roundoff errors would be larger than the theoretical errors (when using $$P_1$$ piecewise affine elements with standard double precision arithmetic)?

Just to fix the framework, I have some particular polygon, meshed with equal triangles and I use $$P_1$$ Lagrange (piecewise affine) elements. Therefore, the stiffness and mass matrices $$K$$ and $$M$$ can be computed explicitly. Zero Dirichlet boundary conditions are applied and the nodes are eliminated from the system by considering submatrices of $$K$$ and $$M$$. When implementing the problem, roundoff errors appear, leading to perturbations $$\Delta K$$ and $$\Delta M$$ of these matrices. In particular, $$M$$ is symmetric, positive definite and $$K$$ is symmetric positive definite (after removal of lines and columns of Dirichlet nodes).

I am interested in solving the generalized eigenvalue problem $$Ku = \lambda Mu.$$ The associated eigenvectors form a basis of $$\Bbb{R}^N$$, orthonormal w.r.t. the scalar product $$x^TMy$$.

Consider now a solution $$\overline u$$ of $$(K+\Delta K)\overline u = \overline \lambda (M+\Delta M)\overline u$$. Write $$\overline u= \sum_{i=1}^N c_i u_i$$ in the basis $$(u_i)$$, and suppose that $$\overline u^T M \overline u=1$$ leading to $$\sum_{i=1}^N c_i^2 = 1$$. After a few computations we have $$\sum_{i=1}^n (\lambda_i-\overline \lambda) c_i M u_i = (-\Delta K+\overline \lambda \Delta M) \overline u$$ Multiply with $$\text{sign}(\lambda_i-\overline \lambda)c_j u_j$$ and sum for $$j=1,...,n$$ to get $$\sum_{i=1}^n c_i^2|\lambda_i-\overline \lambda| = \overline u_m^T(-\Delta K+\overline \lambda \Delta M)\overline u$$ where $$\overline u_m = \sum_{i=1}^N \pm c_i u_i$$ leading to $$\overline u_m^T M \overline u_m=1$$.

Now suppose $$\Delta K$$ and $$\Delta M$$ are machine erorrs for $$K$$ and $$M$$, in particular, componentwise inequalities of the form $$|\Delta K|\leq \varepsilon |K| , |\Delta M|\leq \varepsilon |M|$$ hold (in particular, zero elements in K,M lead to zero elements in $$\Delta K, \Delta M$$)

In the following paper the authors describe how to estimate the eigenvalues of $$K$$ and $$M$$ (maybe other references exist, I found this well explained and nice to read, with quantifiable constants).

In particular, given the mesh size $$h$$, eigenvalues of $$M$$ are all of size $$O(h^2)$$, implying that $$\|\overline u\|_2 = O(h^{-1})$$. Suppose that $$\overline \lambda$$ approximates the smallest eigenvalue $$\lambda_1$$. Then it is straightforward to see (using the fact that the lines of $$K,M$$ have a small number of non-zero elements) that $$|\overline u_m^T \Delta M u_m| \leq \varepsilon O(1), |\overline u_m^T \Delta K u_m| \leq \varepsilon O(h^{-2})$$

Therefore $$\sum_{i=1}^n c_i^2|\lambda_i-\overline \lambda| \leq \varepsilon O(h^{-2})$$ This immediately shows that $$\min |\lambda_i-\overline \lambda| \leq \varepsilon O(h^{-2})$$ Moreover, for $$j>1$$ if $$v_1, \overline v$$ are the $$P_1$$ FEM functions given by coefficients $$u_1$$ and $$\overline u$$ then $$\|v_1-\overline v\|_{L^2}^2 = (\overline u-u_1)^T M (\overline u-u_1) = (c_1-1)^2+\sum_{j>1} c_j^2.$$ $$\|\nabla v_1- \nabla \overline v\|_{L^2}^2 = (\overline u-u_1)^T K (\overline u-u_1) = \lambda_1(c_1-1)^2+\sum_{j>1} \lambda_j c_j^2.$$ The previous estimate also shows that these two quantities are of order $$\varepsilon O(h^{-2})$$.

These error estimates seem too pessimistic. For $$h<1e-4$$ machine precision error would dominate when using standard double precision. Also I would expect the $$L^2$$ norm to be better approximated than the $$H^1$$ one. Is there any way of improving them?

• Have you tried running this problem on smaller and smaller meshes until you hit a problem? You may develop a hint to the theoretical answer you seek. Oct 20 at 14:09
• @BillBarth: Thank you for this suggestion. However, I don't think I can get to a very small $h$ without high performance computing. I will update my questions once my experiments are done. Oct 20 at 23:05
• I only meant on the smallest you could manage. You only need three data points or so in $h$ to get a feel for the first digit and size of the exponent, unless you’re already near machine zero. Oct 20 at 23:11
• @CarlChristian: My goal is to understand roundoff errors. The convergence of the discrete eigenvalues towards the continuous ones is of order $O(h^2)$. However, as $h$ goes to zero, the size of the matrices goes to $\infty$. The computations in the question show that at fixed $h$ in presence of errors of size $\varepsilon$ in the coeffs of the matrices one could guarantee a precision of $\varepsilon O(h^{-2})$ between the eigenvalues of the exact discrete system and the perturbed one. Oct 21 at 20:07
• It is possible that the following this answer might be relevant to your problem. The technique is quite general and should apply to your situation, provided that your functions are sufficiently smooth. Otherwise, the necessary asymptotic error expansions might not exists. It will not allow you estimate the rounding error, but you will be able to say when it is irrelevant compared with the discretization error and you will be able to estimate the discretization error reliably without knowing the exact solution. Oct 21 at 21:11

• Yes, you are right. I wrote my comment before reading the paper thoroughly. Reading it and the references therein helped me understand a bit more clearly the problem. Applying results in Section 3 of the paper lead to estimates of the same order. Repeating point by point the proof in this case gives $\min |\mu-\lambda_i| \leq \|M^{-1/2}EM^{1/2}\| = O(h^{-2})\|E\|$. Oct 21 at 20:13