24. Finite element error analysis#
Let \(u\) be the solution of the variational problem, and \(u_h\) its Galerkin approximation in the finite element sub-space \(V_h\). Cea’s Lemma bounds the finite element error \(u - u_h\) by the best approximation error
The constant factor \(C\) is the ratio of the continuity bound and the coercivity bound of the bilinear form \(A(.,.)\).
Provided that the solution \(u\) is sufficiently smooth, we can take the finite element interpolant to bound the best approximation error:
In the following, we will bound the interpolation error.
Lemma: Let \(\widehat T\) and \(T\) be \(d\)-dimensional domains related by the invertible affine linear transformation \(\Phi_T : \widehat T \rightarrow T \)
\[ \Phi_T (x) = a + B x, \]where \(a \in {\mathbb R}^d\) and \(B\) is a regular matrix in \({\mathbb R}^{d\times d}\). Then there holds:
\[ \| u \|_{L_2(T)} = (\det B)^{1/2} \, \| u \circ \Phi_T \|_{L_2(\widehat T)} \]\[\begin{align*} \frac{\partial}{\partial x_{i_m}} \ldots \frac{\partial}{\partial x_{i_1}} \left( u \circ \Phi_T \right) = \sum_{j_m = 1}^d \ldots \sum_{j_1=1}^d \left( \frac{\partial}{\partial x_{j_m}} \ldots \frac{\partial}{\partial x_{j_1}} u \right) \circ \Phi_T \; \; B_{j_m,i_m} \ldots B_{j_1, i_1} \end{align*}\]\[\begin{equation*} | u \circ \Phi |_{H^m(\widehat T)} \preceq (\det B)^{-1/2} \| B \|^m \, | u |_{H^m(T)} \end{equation*}\]
Proof: Transformation of integrals, chain rule. \(\Box\)
We define the element-size \(h_T\) as the diameter (what is \(\max_{x,y \in T} \| x - y \|\)) of the element \(T\):
A triangulation is called shape regular, if all its elements fulfill
with a good constant \(\simeq 1\). If one studies convergence, one considers families of triangulations with decreasing element sizes \(h_T\). In that case, the family of triangulations is called shape regular, if there is a common constant \(C\) such that all elements of all triangulations fulfill \(| T | \geq C h_T^d\).
Lemma: Let \(\Phi_T = a + B x\) be the mapping from the reference triangle to the triangle \(T\). Let \(|T| \succeq h_T^d\). Then there holds
The following lemma is the basis for the error estimate. This lemma is the main application for the Bramble Hilbert lemma. Sometimes, it is called the Bramble Hilbert lemma itself:
Bramble-Hilbert Lemma:
Let \((T, V_T, \Psi_T)\) be a Lagrange finite element such that the element space \(V_T\) contains polynomials up to order \(P^k\). Then there holds
\[ \| v - I_T v \|_{H^1} \leq C | v |_{H^m} \qquad \forall \, v \in H^m(T) \]for all \(m > d/2\), \(m \geq 1\), and \(m \leq k+1\).
Proof: First, we prove that \(id - I_T\) is a bounded operator from \(H^m\) to \(H^1\):
The last step used that for \(H^m\), with \(m > d/2\), point evaluation is continuous. Now, let \(v \in P^k(T)\). Since \(P^k \subset V_T\), and \(I_T\) is a projection on \(V_T\), there holds \(v - I_T v = 0\). The (first version of the) Bramble Hilbert Lemma applied for \(U = H^1\) and \(L = id - I_T\) proves the result. \(\Box\)
To bound the finite element interpolation error, we will transform functions from the elements \(T\) to the reference element \(\widehat T\).
Theorem: Let \({\cal T}\) be a shape regular triangulation of \(\Omega\). Let \(V_{\cal T}\) be a \(C^0\)-regular finite element space such that all local spaces contain \(P^1\). Then there holds
Proof: We prove the \(H^1\) estimate, the \(L_2\) one follows the same lines. The interpolation error on each element is transformed to the interpolation error on one reference element:
On the reference element \(\widehat T\) we apply the Bramble-Hilbert lemma. Then, we transform back to the individual elements:
\(\Box\)
A triangulation is called \({\bf quasi-uniform}\), if all elements are essentially of the same size, i.e., there exists one global \(h\) such that
On a quasi-uniform mesh, there hold the interpolation error estimates
We are interested in the rate of the error in terms of the mesh-size \(h\).
Theorem: Finite element error estimate
Assume that
the solution \(u\) of the weak bvp is in \(H^2\),
the triangulation \({\cal T}\) is quasi-uniform of mesh-size \(h\),
the element spaces contain \(P^1\).
Then, the finite element error is bounded by
\[ \| u - u_h \|_{H^1} \preceq h \, | u |_{H^2} \]
Proof: the result follows directly from Cea’s lemma and estimates for the interpolation error. \(\Box\)
24.1. Error estimates in \(L_2\)-norm#
The above theorem bounds the error in the \(L_2\)-norm of the function, and the \(L_2\)-norm of the derivatives with the same rate in terms of \(h\). This is obtained by the natural norm of the variational formulation.
The interpolation error suggests a faster convergence in the weaker norm \(L_2\). Under certain circumstances, the finite element error measured in \(L_2\) also decays faster. The considered variational problem is
We define the dual problem as
In the case of a symmetric bilinear form, the primal and the dual problem coincide.
Aubin-Nitsche Trick:
Assume that
the dual weak bvp is \(H^2\) regular
the triangulation \({\cal T}\) is quasi-uniform of mesh-size \(h\),
the element spaces contain \(P^1\). Then, there holds the \(L_2\)-error estimate
\[ \| u - u_h \|_{L_2} \preceq h^2 \, | u |_{H^2} \]
Proof: Solve the dual problem with the error \(u - u_h\) as right hand side:
Since the dual problem is \(H^2\) regular, there holds \(w\in H^2\), and \(\| w \|_{H^2} \preceq \| u-u_h \|_{L_2}\). Choose the test function \(v := u - u_h\) to obtain the squared norm
Using the Galerkin orthogonality \(A(u-u_h, v_h) = 0\) for all \(v_h \in V_h\), we can insert \(I_{\cal T} w\):
Next we use continuity of \(A(.,.)\) and the interpolation error estimates:
From \(H^2\) regularity:
and, after dividing one factor
\(\Box\)
24.2. Approximation of Dirichlet boundary conditions#
Till now, we have neglected Dirichlet boundary conditions. In this case, the continuous problem is
where
The finite element problem is
where
The definition of \(V_{hD}\) coincides with \( \{ v_h \in V_h : v_h(x_i) = u_D(x_i) \; \; \forall \mbox{ vertices } x_i \mbox{ on } \Gamma_D \}. \)
There holds \(V_{h0} \subset V_0\), but, in general, there does not hold \(V_{hD} \subset V_D\).
Error estimate for Dirichlet boundary conditions Assume that
\(A(.,.)\) is coercive on \(V_{h0}\):
\[ A(v_h, v_h) \geq \alpha_1 \, \| v_h \|_V^2 \qquad \forall \, v_h \in V_{h0} \]
\(A(.,.)\) is continuous on \(V\):
\[ A(u, v) \leq \alpha_2 \, \| u \|_V \| v \|_V \qquad \forall \, u, v \in V \]Then there holds the finite element error estimate
\[ \| u - u_h \|_{H^1} \preceq h | u |_{H^2} \]
Proof: To make use of the coercivity of \(A(.,.)\), we need an element in \(V_{h0}\). There holds Galerkin orthogonality \(A(u-u_h, v_h) = 0 \; \forall \, v_h \in V_{h0}\):
Next, we apply \(a b \leq \frac{1}{2} a^2 + \frac{1}{2} b^2\) for \(a = \frac{2 \alpha_2}{\alpha_1} \| u - I_h u \|_V\) and \(b = \| u - u_h \|_V\):
Moving the term \(\frac{1}{2} \| u - u_h \|\) to the left, we obtain
\(\Box\)
24.3. High order elements#
One can obtain faster convergence, if the solution is smooth, and elements of higher order are used:
Theorem: Assume that
the solution is smooth: \(u \in H^m\) for \(m \geq 2\)
all element spaces \(V_T\) contain polynomials \(P^p\) for \(p \geq 1\)
the mesh is quasi-uniform
Then there holds
\[ h^{-1} \| u - I_h u \|_{L_2} + \| u - I_h u \|_{H^1} \preceq h^{\min\{ m-1,p\} } \| u \|_{H^m} \]
The proof is analogous to the case \(m=2\) and \(p=1\). The constants in the estimates depend on the Sobolev index \(m\) and on the polynomial order \(p\). Nodal interpolation is instable (i.e., the constant grow with \(p\)) for increasing order \(p\). There exist better choices to bound the best approximation error.
24.4. Graded meshes around vertex singularities#
On non-convex domains, the solution is in general not in \(H^2\), but in some weighted Sobolev space. The information of the weight can be used to construct proper locally refined meshes.
On a sector domain with a non-convex corner of angle \(\omega > \pi\), the solution is bounded in the weighted Sobolev norm
with \(\beta = \frac{\pi}{\omega} \in (\tfrac{1}{2}, 1)\). One may choose a mesh such that
where \(r_T\) is the distance of the center of the element to the singular corner, and \(\underline h \in {\mathbb R}^+\) is a global mesh size parameter.
We bound the interpolation error:
The number of elements in the domain can be roughly estimated by the integral over the density of elements. The density is number of elements per unit volume, i.e., the inverse of the area of the element:
In two dimensions, and \(\beta \in (0,1)\), the integral is finite.
Combining the two estimates, one obtains a relation between the error and the number of elements:
This is the same order of convergence as in the \(H^2\) regular case !
Building the meshes according to the weighted regularity requires a lot of knowledge about the solution. Adaptive refinement algorithms recover quasi-optimal meshes automatically.
24.5. Geometric mesh refinement#
Together with high order finite element spaces, geometric mesh refinement towards singular points leads to an optimal refinement strategy. Only the elements sitting next to the vertex are sub-divided. Convergence is obtained by simultaneously increasing the polynomial order \(p\). Together it is called \(hp\)-refinement.
from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
shape = MoveTo(0,0).Rotate(30).Line(1).Rotate(90).Arc(1, 305).Close().Face()
shape.vertices.Nearest((0,0)).hpref=1
mesh = Mesh(OCCGeometry(shape,dim=2).GenerateMesh(maxh=0.5)).Curve(3)
mesh.RefineHP(levels=3, factor=0.25)
Draw(mesh);
Close to edges in 3D, the solution is smooth along the edge, but derivatives normal to the edge are singular. This requires anisotropic meshes, which have small element size across the edge, but large element size along the edge. A careful analysis combines small sizes and large derivatives in different directions. By this refinement, meshes with different element shapes are obtained:
cyl1 = Cylinder(Axes(), r=2, h=1)
cyl2 = Cylinder(Axes(), r=1, h=2)
cyl1.faces.name="cyl1"
cyl2.faces.name="cyl2"
shape = cyl1+cyl2
(shape.faces["cyl1"].edges*shape.faces["cyl2"].edges).hpref=1
mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.5)).Curve(3)
mesh.RefineHP(levels=2, factor=0.25)
Draw (mesh, order=3);