30. hp - Finite Elements

Let \(V_{hp}\) be a \(p\)-th order finite element sub-space of \(H^1\). By scaling and Bramble-Hilbert technique one obtains the best-approxiamtion error estimate

\[ \inf_{v_{hp} \in V_{hp}} \| u - v_{hp} \|_{H^1} \leq c h^{m-1} \| u \|_{H^m} \]

for \(m \leq p+1\). The constant \(c\) depends on the order \(p\). If \(m\) is fixed, we do obtain reduction of the approximation error as we increase \(p\). Next we develop methods to obtain so called \(p\)-version error estimates

\[ \inf_{v_{hp} \in V_{hp}} \| u - v_{hp} \|_{H^1} \leq c \left( \tfrac{h}{p} \right) ^{m-1} \| u \|_{H^m}, \]

where \(c\) is independent of \(h\) and \(p\). This estimate proves also convergence of the \(p\)-version finite element method: One may fix the mesh, and increase the order \(p\).

A detailed analysis of local \(H^m\) norms allows an optimal balance of mesh-size \(h\) and polynomial order \(p\). This \(hp\)-version leads to exponential convergence

\[ \inf_{v_{hp} \in V_{hp}} \| u - v_{hp} \|_{H^1} \leq c e^{ -N^\alpha }, \]

where \(N\) is the number of unknowns.

We will prove the \(p\)-version estimate, but not the \(hp\)-result.

30.1. Legendre Polynomials

Orthogonal polynomials are important to construct stable basis functions for the \(p\)-FEM, as well as for error estimates.

Let \(\Pi_n\) denote the space of polynomials up to order \(n\). We write \(\pi_n\) for a generic polynomial in \(\Pi_n\), with a different value any time it appears.

Definition of Legendre polynomials via Rodrigues’ formula:

\[ P_n(x) := \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2-1)^n. \]

It is a polynomial of degree \(n\). The first few Legendre polynomials are

\[\begin{align*} P_0(x) & = 1 \\ P_1(x) & = x \\ P_2(x) & = \tfrac{3}{2} x^2 - \tfrac{1}{2} \end{align*}\]

\(P_n\) is even if \(n\) is even, and \(P_n\) is odd if \(n\) is odd. Since \((x^2-1)^n = x^{2n} - n x^{2n-2} + \pi_{2n-4}\) (with proper modification for small \(n\)) we have

\[ % \label{equ_leadingcoef} P_n(x) = \frac{1}{2^n n!} \frac{ (2n)! }{ n! } x^n - \frac{n}{2^n n!} \frac{ (2n-2)!} {(n-2)!} x^{n-2} + \pi_{n-4} \]

Lemma: Orthogonality

There holds

\[ \int_{-1}^1 P_n(x) P_m(x) \, dx = \frac{2}{2n+1} \delta_{n,m}. \]

Proof: W.l.o.g. let \(n \leq m\). Multiple integration by parts gives

\[\begin{align*} 2^{n+m} n! m! \int_{-1}^1 P_n(x) P_m(x) \, dx &= \int_{-1}^1 \frac{d^n}{dx^n} (x^2-1)^n \frac{d^m}{dx^m} (x^2-1)^n \, dx \\ & = \int_{-1}^1 \frac{d^{n+1}}{dx^{n+1}} (x^2-1)^n \frac{d^{m-1}}{dx^{m-1}} (x^2-1)^m + \left[ \frac{d^n}{dx^{n}} (x^2-1)^n \underbrace{ \frac{d^{m-1}}{dx^{m-1}} (x^2-1)^m}_{= 0 \text{ for } x \in \{ -1, 1\}} \right]_{-1}^1 \\ & = \cdots \\ & = \int_{-1}^1 \frac{d^{n+m}}{dx^{n+m}} (x^2-1)^n (x^2-1)^m \, dx \end{align*}\]

For \(n < m\), the first factor of the integrand vanishes, and we have orthogonality. For \(n = m\) this equals

\[\begin{eqnarray*} 2^{2n} (n!)^2 \| P_n \|_{L_2}^2 &= &\int_{-1}^1 (2n)! (x^2-1)^n \, dx = (2n)! \int_{-1}^1 (x-1)^n (x+1)^n \\ & = & -(2n)! \int_{-1}^1 \frac{n}{n+1} (x-1)^{n+1} (x+1)^{n-1} \\ & = & (2n)! \int_{-1}^1 \frac{ n (n-1) }{ (n+1) (n+2) } (x-1)^{n+2}(x+1)^{n-2} = ... \\ & = & (2n)! \frac{n!} {2n (2n-1) \cdots (n+1)} \int_{-1}^1 (x-1)^{2 n} \, dx = (n!)^2 \frac{1}{2n+1} \, 2^{2n+1}, \end{eqnarray*}\]

which proves the scaling. \(\Box\)

Next we prove the 3-term recurrency, which can be used for efficient evaluation.

Lemma: Three term recurrence:

There holds

(30.1)\[\begin{equation} % \label{equ_threeterm} (n+1) P_{n+1} (x) = (2n+1) x P_n(x) - n P_{n-1}(x). \end{equation}\]

Proof: Set \(r(x) = (n+1) P_{n+1} (x) - (2n+1) x P_n(x) + n P_{n-1}(x)\). Using (\ref{equ_leadingcoef}), we see that the leading coefficients cancel, and thus \(r \in \Pi_{n-2}\). From Lemma~\ref{lemma_ortho} we get for any \(q \in \Pi_{n-2}\)

\[ \int_{-1}^1 r(x) q(x) \, dx = (n+1) \int_{-1}^1 P_{n+1} \, q - (2n+1) \int_{-1}^1 P_n \underbrace{x q}_{\in \Pi_{n-1}} + n \int_{-1}^1 P_{n-1} \, q = 0, \]

and thus \(r = 0\). \(\Box\)

Lemma: Sturm-Liouville

Legendre polynomials satisfy the Sturm-Liouville differential equation

\[ \frac{d}{dx} \big[ (x^2-1) \frac{d}{dx} P_n(x) \big] = n (n+1) P_n(x) \]

Proof: Both sides are in \(\Pi_n\). We compare leading coefficients, for this set \(P_n = a_n x^n + \pi_{n-2}\) (with \(a_n = \frac{1}{2^n n!} \frac{ (2n)! }{ n! }\)).

\[\begin{eqnarray*} lhs & = & \frac{d}{dx} \left[ (x^2-1) \frac{d}{dx} \left( a_n x^n + \pi_{n-2} \right) \right] \\ & = & \frac{d}{dx} \left[ (x^2-1) (a_n n x^{n-1} + \pi_{n-3}) \right] \\ & = & \frac{d}{dx} \left[ a_n n x^{n+1} + \pi_{n-1} \right] \\ & = & n (n+1) a_n x^n + \pi_{n-2}, \end{eqnarray*}\]

and we get the same leading coefficient for rhs. Furthermore, for \(q \in \Pi_{n-1}\) there holds

\[\begin{eqnarray*} \int_{-1}^1 lhs \, q & = &-\int_{-1}^1 (x^2-1) P_n^\prime q^\prime \, dx + \underbrace{ \left[ (x^2-1) P_n^\prime q \right]_{-1}^1 }_{= 0} \\ & = & \int_{-1}^1 P_n \underbrace{ \left( (x^2-1) q^\prime\right)^\prime}_{\in \Pi_{n-1}} \, dx - \left[ P_n (x^2-1) q^\prime \right]_{-1}^1 = 0, \end{eqnarray*}\]

and the same for the rhs. Thus the identity is proven. \end{proof}

Lemma~\ref{lemma_sturmliouville} implies that the Legendre polynomials are also orthogonal w.r.t. \((u^\prime, v^\prime)_{L_2, 1-x^2}\), i.e. $\( \int_{-1}^1 (1-x^2) P_n^\prime P_m^\prime = n(n+1) \| P_n \|_{L_2}^2 \delta_{n,m} \)$

30.2. Error estimate of the \(L_2\) projection

Since polynomials are dense in \(L_2(-1,1)\), we get

\[ u = \sum_{n=0}^\infty a_n P_n \]

with the generalized Fourier coefficients

\[ a_n = \frac{ (u, P_n)_{L_2} } { \| P_n \|^2_{L_2} }, \]

and

\[ \| u \|_{L_2}^2 = \sum_{n=0}^\infty a_n^2 \| P_n \|^2. \]

Let \(P_{L_2}^{\Pi_p}\) denote the \(L_2\)-projection onto \(\Pi_p\). There holds

\[ P_{L_2}^{\Pi_p} u = \sum_{n=0}^p a_n P_n \]

The projection error is

\[ \| u - P_{L_2}^{\Pi_p} u \|_{L_2}^2 = \sum_{n=p+1}^\infty a_n^2 \| P_n \|^2 \]

Lemma Projection error:

The \(L_2\)-projection error satisfies

(30.2)\[\begin{equation} \| u - P_{L_2}^{\Pi_p} u \|_{L_2(-1,1)} \leq \frac{1}{\sqrt{(p+1)(p+2)}} \, | u |_{H^1(-1,1)} \end{equation}\]

Proof: Since \(P_n\) are orthogonal also w.r.t. \((u^\prime, v^\prime)_{L_2, 1-x^2}\), there holds

\[ \| u^\prime \|_{L_2, 1-x^2}^2 = \sum_{n \in {\mathbb N}} a_n^2 \, \| P_n^\prime \|_{1-x^2}^2, \]

provided that \(u\) is in \(H^1\). The projection error satisfies

\[\begin{align*} \| u - P_{L_2}^{\Pi_p} u \|^2_{L_2} & = \sum_{n > p} a_n^2 \| P_n \|^2_{L_2} = \sum_{n > p} a_n^2 \frac{1}{n (n+1)} \| P_n^\prime \|^2_{1-x^2} \\ & \leq \frac{1}{(p+1)(p+2)} \sum_{n > p} a_n^2 \|P_n^\prime \|^2_{1-x^2} \leq \frac{1}{(p+1)(p+2)} \sum_{n \in {\mathbb N}} a_n^2 \|P_n^\prime \|^2_{1-x^2} \\ & = \frac{1}{(p+1)(p+2)} \| u^\prime \|_{1-x^2}^2 \\ \end{align*}\]

Finally, the result follows from

\[ \int_{-1}^1 (1-x^2) (u^\prime)^2 \, dx \leq \int_{-1}^1 (u^\prime)^2 \, dx. \]

\(\Box\)

Similar as in Lemma~\ref{lemma_sturmliouville} on shows also

\[ \frac{d^m}{dx^m} \big[ (x^2-1)^m \frac{d^m}{dx^m} P_n(x) \big] = (n+m)(n+m-1) \ldots (n-m+1) P_n(x) \]

for \(m \leq n\), and, as in Lemma~\ref{lemma_l2est}

\[ \| u - P_{L_2}^{\Pi_p} u \|_{L_2} \leq \sqrt{ \frac{ (p-m+1)! }{ (p+m+1)! } } | u |_{H^m}. \]

30.3. Orthogonal polynomials on triangles

Orthogonal polynomials on tensor product elements are simply constructed by tensorization. Orthogonal polynomials on simplicial elements are more advanced. They are based on Jacobi polynomials:

For \(\alpha, \beta > -1\), Jacobi polynomials are defined by

\[ P_n^{(\alpha, \beta)}(x) := \frac{(-1)^n}{2^n n!} \frac{1}{w(x)} \frac{d^n}{dx^n} \big( w(x) (1-x^2)^n \big) \]

with the weight function

\[ w(x) = (1-x)^\alpha (1+x)^\beta. \]

Jacobi polynomials are orthogonal w.r.t. the weighted inner product

\[ \int_{-1}^1 w(x) P_n^{(\alpha,\beta)}(x) P_m^{(\alpha, \beta)}(x) dx = \delta_{n,m} \frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1} \frac{ \Gamma(n+\alpha+1) \, \Gamma (n+\beta+1) } { n! \Gamma (n+\alpha+\beta+1)}. \]

Note that \(P^{(0,0)}_n = P_n\).

Define the unit-triangle \(T\) with vertices \((-1,0)\), \((1,0)\) and \((0,1)\).

Lemma: Dubiner basis

The functions

\[ \varphi_{i,j} (x,y) := P_i \big( \tfrac{x}{1-y} \big) (1-y)^i P_j^{(2i+1,0)} (2y-1) \qquad i+j \leq p \]

form an \(L_2(T)\)-orthogonal basis for \(\Pi_p(T)\).

Proof: Note that \(\varphi_{i,j} \in \Pi_{i+j}\). Substitution \(\xi = \frac{x}{1-y}\) leads to

\[\begin{eqnarray*} \int_T \varphi_{ij}(x,y) \varphi_{kl} (x,y) \, d(x,y) &= & \int_0^1 \int_{-1+y}^{1-y} P_i \big( \tfrac{x}{1-y} \big) (1-y)^i P_j^{(2i+1,0)} (2y-1) P_k \big( \tfrac{x}{1-y} \big) (1-y)^k P_l^{(2k+1,0)} (2y-1) \, dx dy \\ & = & \int_0^1 \int_{-1}^1 P_i(\xi) P_k(\xi) (1-y)^{i+k+1} P_j^{(2i+1,0)} (2y-1) P_l^{(2k+1,0)}(2y-1) \, d\xi dy \\ & = & \delta_{i,k} \| P_i \|_{L_2}^2 \int_0^1 (1-y)^{2i+1} P_j^{(2i+1,0)}(2y-1) P_l^{(2i+1,0)}(2y-1) \, dy \\ & = & C_{ij} \delta_{i,k} \delta_{j,l} \end{eqnarray*}\]

\(\Box\)

30.4. Projection based interpolation

By means of the orthogonal polynomials one shows approximation error estimates of the form $\( \inf_{q \in \Pi_p(T)} \| u - q \|_{H^k(T)} \leq c p^{k-m} | u |_{H^m(T)} \qquad m \geq k, \)\( with \)c \neq c(p)\(, easily in 1D and tensor product elements, and also on \)n$-dimensional simplices [Braess+Schwab: Approximation on simplices with respect to weighted Sobolev norms, J. Approximation Theory 103, 329-337 (2000)].

But, an interpolation operator to an \(H^1\)-conforming finite element space has to satisfy continuity constraints across element boundaries. We show that we get the same rate of convergence under these constraints.

30.4.1. The 1D case

Let \(I=(-1,1)\). We define the operator \(I_p : H^1(I) \rightarrow \Pi_p\) such that

\[\begin{align*} I_p u(x) & = u(x) \qquad x \in \{ -1, 1 \} \\ \int_I (I_p u)^\prime q^\prime & = \int_I u^\prime q^\prime \qquad \forall \, q \in \Pi_{p,0}(I), \end{align*}\]

where \(\Pi_{p,0}(D) := \{ q \in \Pi_p(D) : q = 0 \text{ on } \partial D\}\). This procedure is exactly a \(p\)-version Galerkin-method for the Dirichlet problem. Since boundary values are preserved, the interpolation operator produces a globally continuous function. The operator \(I_p\) is a kind of mixture of interpolation and projection, thus the term {\em projection based interpolation} introduced by Demkowicz has been established.

Lemma: Commuting diagram There holds

\[ \Pi_{L_2}^{\Pi_{p-1}} u^\prime = (I_p u)^\prime \]

Proof: The range of both sides is \(\Pi_{p-1}\). We have to show that \((I_p u)^\prime\) is indeed the \(L_2\)-projection of \(u^\prime\), i.e.

\[ \int_I (I_p u)^\prime q = \int_I u^\prime q \qquad \forall \, q \in \Pi_{p-1}. \]

This holds since \(\{ q^\prime : q \in \Pi_{p,0} \} = \{ q \in \Pi_{p-1} : \int q = 0 \}\) and (\ref{equ_projbased2}), and

\[ \int_{-1}^1 (I_p u)^\prime \, 1 = (I_p u)(1) - (I_p u)(-1) = u(1) - u(-1) = \int_{-1}^1 u^\prime \, 1. \]

\(\Box\)

The \(H^1\)-error estimate follows directly from the commuting diagram property:

\[ | u - I_p u |_{H^1(I)} = \| u^\prime - (I_p u)^\prime \|_{L_2} = \| u^\prime - P_{L_2}^{\Pi_{p-1}} u^\prime \|_{L_2} \leq \frac{c}{p^{m-1}} | u^\prime |_{H^{m-1}}. \]

By the Aubin-Nitsche technique one obtains an extra \(p\) for the \(L_2\)-error:

\[ \| u - I_p u \|_{L_2(I)} \preceq \frac{1}{p} | u - I_p u |_{H^1(I)} \leq \frac{1}{p^m} | u |_{H^m(I)} \]

One also gets for \(q \in \Pi_p\)

\[ | u - I_p u |_{H^1} = | u-q - I_p (u-q) |_1 \leq || Id - I_p ||_{H^1\rightarrow H^1} \, | u - q |_{H^1} \]

30.4.2. Projection based interpolation on triangles

We define the operator \(I_p : H^2(T) \rightarrow \Pi_p(T)\) as follows:

(30.3)\[\begin{eqnarray} I_p u(x) & = & u(x) \qquad \forall \text{ vertices } x \\ \int_E \partial_\tau (I_p u) \partial_\tau q & = & \int_E \partial_\tau u \partial_\tau q \qquad \forall \text{ edges } E, \; \forall \, q \in \Pi_{p,0}(E) \\ \int_T \nabla (I_p u) \nabla q & = & \int_T \nabla u \nabla q \qquad \forall \, q \in \Pi_{p,0}(T) \end{eqnarray}\]

Note that \(I_p u\) on the edge \(E\) depends only on \(u|_E\), and thus the interpolant is continuous across neighbouring elements.

Lemma: Let \(v \in C(\partial T)\) such that \(v|_E \in \Pi_p(E)\). Then there exists an extension \(\tilde v \in \Pi_p(T)\) such that \(\tilde v|_{\partial T} = v\) and

\[ | \tilde v |_{H^1} \leq c \, | v |_{H^{1/2}(\partial T)}, \]

where \(c\) is independent of \(p\).

Major steps have been shown in exercises 5.2 and 6.6. Note that the minimal-norm extension \(\tilde v\) is the solution of the Dirichlet problem, i.e. $\( \int_T \nabla \tilde v \nabla w = 0 \qquad \forall \, w \in \Pi_{p,0} \)$

Theorem: error estimate

There holds

\[ | u - I_p u |_{H^1} \preceq \inf_{q \in \Pi_p} | u - q |_{H^1(T)} + \sum_{E \subset \partial T} \frac{1}{\sqrt{p}} \inf_{q \in \Pi_p(E)} | u - q |_{H^1(E)} \preceq \frac{1}{p^{m-1}} | u |_{H^m} \]

for \(u \in H^m, m \geq 2\).

Proof: Let \(u_p\) be the \(|\cdot |_{H^1}\) best approximation to \(u\), i.e.

\[ \int_T \nabla u_p \nabla v = \int_T \nabla u \nabla v \qquad \forall \, v \in \Pi_p, \]

and, for uniqueness, mean values are preserved: \(\int_T u_p = \int_T u\). There holds

\[ | u - u_p |_{H^1} \leq \frac{c}{p^{m-1}} | u |_{H^m}. \]

We apply the triangle inequality:

\[ | u - I_p u |_{H^1} \leq | u - u_p |_{H^1} + | u_p - I_p u |_{H^1} \]

Since

\[ \int_T \nabla u_p \nabla v = \int_T \nabla u \nabla v = \int_T \nabla I_p u \nabla v \qquad \forall \, v \in \Pi_{p,0}(T), \]

we have that

\[ u_p - I_p u \; \bot_{H^1} \; \Pi_{p,0}, \]

i.e. \(u_p - I_p u\) is solution of the Dirichlet problem with boundary values \( (u_p - I_p u)|_{\partial T}\). Lemma~\ref{lemma_polext} implies that

\[ | u_p - I_p u |_{H^1(T)} \preceq | u_p - I_p u |_{H^{1/2}(\partial T)} \]

We insert an \(u\) on the boundary to obtain

\[\begin{eqnarray*} | u - I_p u |_{H^1} & \preceq & | u - u_p |_{H^1(T)} + | u_p - I_p u |_{H^{1/2}(\partial T)} \\ & \leq & | u - u_p |_{H^1(T)} + | u_p - u |_{H^{1/2}(\partial T)} + | u - I_p u |_{H^{1/2}(\partial T)} \\ & \leq & | u - u_p |_{H^1(T)} + \| u - I_p u \|_{L_2(\partial T)}^{1/2} | u - I_p u |_{H^1(\partial T)}^{1/2}. \end{eqnarray*}\]

In the last step we used that \(H^{1/2}(\partial T) = [L_2, H^1]_{1/2}\) (i.e. the interpolation space). Next, we observe that \(I_p\) restricted to one edge \(E\) is exactly the 1D operator. Using Aubin-Nitsche we get

\[\begin{eqnarray*} | u - I_p u |_{H^1} & \preceq & | u - u_p |_{H^1(T)} + p^{-1/2} \| u - I_p u\|_{H^1(\partial T)} \\ & \preceq & | u - u_p |_{H^1(T)} + \sum_E p^{1-m} \, | u |_{H^{m-1/2}(E)} \\ & \preceq & p^{1-m} | u |_{H^m(T)} \end{eqnarray*}\]

In the last step we used the trace theorem. \(\Box\)