How can I derive the a priori error estimate for a symmetric bilinear form using lagrange finite elements?

Question

Suppose that I'm solving the poisson equation by the finite element method by lagrange elements. I know that the error can be measured in a variety of ways, depending on which norm you choose. For a given lagrange element order p and a symmetric bilinear form (e.g. galerkin's finite element method), the errors in the $L^2$ and $H^1$ norms are

$$||u-u_h||_{L^2}\le Ch^{p+1} |u|_{H^{p+1}}$$ and $$||u-u_h||_{H^1}\le Ch^p|u|_{H^{p+1}}$$

where $|u|_{H^{p+1}}$ is the seminorm $$|u|^2_{H^{p+1}} = \sum_{i+j=p+1}\int_\Omega |\frac{\partial^{p+1}u}{\partial x^i\partial y^j}|^2 $$

I before I begin delving into reading the full details of the proof of this result (e.g. Brenner & Scott (1994) The Mathematical Theory of Finite Element Methods, Chapter 4), I want to get a "big picture" sketch of the derivation of the bounds. What are the key ideas that go into the proof and how do they fit together to bound the error?

My hope is to gain insight into the overall strategy of deriving error bounds on other (perhaps more exotic) finite elements by first understanding how to derive them for a well-studied case such as this.

Answer 1

To extend VorKir's basically correct answer, to derive a priori error estimate for finite element approximations, you need to stack the following three Lego blocks:

1. The well-posedness of the weak formulation of the PDE.
2. The (global) conformity of the Galerkin approximation.
3. The (local) approximation properties of your finite element space.

(Only the last one is specific to finite elements.) Of course, you need to stack them correctly, so let's look at these in a bit more detail. I'll state them abstractly since you're interested in the blueprint and not the specifics.

Let's assume your problem in weak formulation is to find $u\in V$ such that $$a(u,v) = (f,v) \qquad\text{for all }v\in V \label{weak}\tag{1}$$ for some (not necessarily symmetric) bilinear form $a:V\times V\to\mathbb{R}$. Well-posedness here means that there exist constants $c,C>0$ such that for all $u,v\in V$, $$\begin{aligned} c\|u\|_V^2 &\leq a(u,u), \\ a(u,v)&\leq C\|u\|_V\|v\|_V.\end{aligned}$$ The first inequality is called coercivity, the second continuity. The space $V$ in which you can prove both is the natural space for the problem and determines the rest of the analysis. Usually $V$ is a Sobolev space such as $H^1(\Omega)$ (which is the standard space for elliptic second order equations). From this, existence, uniqueness and stability of the solution follows using the Lax-Milgram Lemma, which is basically a nonsymmetric variant of the Riesz representation theorem. (If your test function needs to be taken in a different space $W\neq V$ (because you can't show continuity otherwise), the coercivity is replaced by the so-called inf-sup inequality, and the Lax-Milgram Lemma by the Banach-Nečas-Babuška Theorem.)

Now you pick a finite-dimensional subspace $V_h\subset V$ (this is the key property!) and look for $u_h\in V_h$ such that $$a(u_h,v_h) = (f,v_h) \qquad\text{for all }v_h\in V_h \label{galerkin}\tag{2}.$$ Since $v_h\in V_h\subset V$, you can insert $v_h$ in \eqref{weak} as well and subtract the two to obtain the Galerkin orthogonality $$a(u-u_h,v_h) = 0 \qquad\text{for any }v_h\in V_h.$$ Hence, for arbitrary $v_h\in V_h$, you have $v_h-u_h\in V_h$ and therefore $a(u-u_h,v_h-u_h)=0$. Using coercivity and continuity of the bilinear form, you obtain $$ \begin{aligned} c \|{u-u_h}\|_V^2& \leq a(u-u_h,u-u_h)\\ &=a(u-u_h,u-v_h) + a(u-u_h,v_h-u_h)\\ &\leq C \|{u-u_h}\|_V\|{u-v_h}\|_V. \end{aligned} $$ Dividing by $\|{u-u_h}\|_V$, rearranging, and taking the infimum over all $v_h\in V_h$ yields $$ \|{u-u_h}\|_V \leq \frac{C}{c} \inf_{v_h\in V_h}\|{u-v_h}\|_V.\label{cea}\tag{3}$$ This is Céa's lemma, which says that the error of any (conforming) Galerkin approach is determined by the approximation error of the exact solution in $V_h$. (If the discrete space is not conforming, i.e., $V_h\not\subset V$, you get extra terms you need to estimate and end up with Strang's Lemmas instead.)

The next step is fairly straightforward: By Céa's Lemma, the discretization error is bounded by the best approximation error in $V$, which in turn can be bounded by any approximation error -- for example, the interpolation error in $V$, which as pointed out above is usually a Sobolev space. The interpolation error of polynomials in Sobolev spaces can be estimated by the Bramble-Hilbert lemma, which yields the error estimate $$\|u-\mathcal{I}_h u\|_{H^m} \leq C h^{k-m} \|u\|_{H^k}\label{poly}\tag{4}$$ provided $V_h$ contains piecewise polynomials of degree $k-1$ (and of course $u\in H^k$, which is an assumption on the solution of \eqref{weak}). (Basically, this follows from the fact that the interpolation error is bounded on the reference element, combined with a simple scaling using transformation of integrals. Summing over all elements gives a global interpolation estimate.) Note that the choice of $m$ in \eqref{poly} is determined by the condition that $H^m\hookrightarrow V$ so that you can chain \eqref{cea} and \eqref{poly} to obtain the standard a priori error estimate in the $V$ norm (I'm being free here with the constants, which can change from one occurrence to the next): $$\label{apriori}\tag{5} \|{u-u_h}\|_V \leq C \inf_{v_h\in V_h}\|{u-v_h}\|_V \leq C \|u-\mathcal{I}_h u\|_{V} \leq C \|u-\mathcal{I}_h u\|_{H^m} \leq C h^{k-m} \|u\|_{H^k}.$$

To get the $L^2$ estimate from this requires the Nitsche trick. Assume you have a Hilbert space $H$ with $V\hookrightarrow H$ (such as $H=L^2(\Omega)$ for $V=H^m(\Omega)$). Using the fact that you can identify the Hilbert space $H$ with its own dual space and thus write the $H$ norm as an operator norm, you can derive the inequality $$\label{aubin}\tag{6} \|u-u_h\|_{H} = \sup_{\varphi\in H\setminus\{0\}} \frac{(u-u_h,\varphi)_H}{\|\varphi\|_{H}} \leq C \|u-u_h\|_{V} \sup_{\varphi\in H\setminus\{0\}} \frac{\inf_{w_h\in V_h}\|w_\phi-w_h\|_{V}}{\|\varphi\|_{H}}, $$ where $w_\varphi$ is the solution of the dual (or adjoint) problem $$\label{adjoint}\tag{7} a(v,w_\varphi) = (\varphi,v)_H \qquad\text{for all }v\in V $$ (note that the test function $v$ is now the first argument of $a$), which is assumed to be well-posed, too. This is the Aubin-Nitsche Lemma (see my previous answer for a derivation in a specific case). The idea is now to use the same interpolation error estimate for the dual problem as well to get an additional power of $h$. First, you have as before $$ \inf_{w_h\in V_h}\|{w_\varphi-w_h}\|_V \leq \|{w_\varphi - \mathcal{I}_h w_\varphi}_V \| \leq C h\|w_\varphi\|_{{H}^2} \leq Ch \|{\varphi}\|_{H}. $$ (In the last step, we have used the well-posedness of the dual problem to estimate the $H^2$ norm of the solution $w_\varphi$ by the $H=L^2$ norm of the right-hand side $\varphi$, which requires some assumptions on the PDE.) Combining this with \eqref{aubin} and \eqref{apriori} gives the estimate $$\|u-u_h\|_{H} \leq C h \|u-u_h\|_{V} \leq C h^{k-m+1} \|u\|_{H^k}.$$

---

So, to summarize, if you want to create your own exotic finite element spaces and show error estimates, you need to

1. Find the correct energy space $V$. (For the Poisson equation with homogeneous Dirichlet conditions, $V=H^1_0(\Omega)$.)

2. Choose the discrete approximation space $V_h$ such that discrete functions in $V_h$ have the right global continuity properties to be in the $V$ (or you can estimate the additional terms from Strang's Lemma). (For Poisson's equation, having $v_h$ continuously differentiable on each element and continuous across element boundaries to ensure $v_h\in V$.)

3. Decide on the polynomial degree you want to represent to get the correct power of $h$, and then choose your degrees of freedom to a) have a unique solution of the interpolation problem in each element and b) ensure the global continuity properties from step 2. (For the Poisson equation and the estimate you gave, this leads to piecewise linear polynomials, and the Lagrange elements correspond to choosing the degrees of freedom at the vertices of the element. Taking $m=1$ and $k=2$ in the interpolation error gives the estimate you stated.)

4. The Nitsche trick then gives you a free additional power of $h$ if you can show that the dual problem is well-posed and you only measure the error in a weaker norm. (For the Poisson equation which is symmetric, the dual problem is of course another Poisson equation, and that step is a freebie.)

---

One final remark: In your case, I would recommend Braess's book over Brenner and Scott's; it's just as mathematically rigorous, but it's meant as a textbook instead of a monograph and thus concentrates on the fundamental ideas rather than the more general and technical results:

Dietrich Braess, MR 2322235 Finite elements, ISBN: 978-0-521-70518-9.

^{(I also can't resist mentioning my lecture notes.)}

Answer 2

I am not pretending to give the best point of view, but I think of the following general lines.

To obtain the error estimates in classical cases, you use as ingredients:

1. interpolation error (which is purely the property of the chosen finite element space and deals with interpolation of one class of functions by another) $$ \parallel u - \pi_h u \parallel \leq Ch^k |u|_{k+1} $$

2. Céa's lemma, which says that the approximate solution is the projection in energy norm of the exact solution onto the finite dimensional subspace. $$ \parallel u - u_h \parallel_A = \inf_{v_h \in V_h} \parallel u - v_h \parallel_A $$ where $V_h$ is the finite element space and $\parallel \cdot \parallel_A$ is the energy norm produced by the main symmetric bilinear form of the problem.

3. Equivalence of the norms in the energy space and appropriate Sobolev space (for the simplest case of Dirichlet problem, it is $H^1$).

This way one can get an estimate in $H^1$. For $L_2$ Nitsche's trick is usually used.