17  FEM: Error Estimates

Lecture 17

Note

I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.

Warning

These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Tip

This chapter is mostly based on

• Burden, Faires, and Burden (2015) Chapter 11,
  
• Driscoll and Braun (2018) Chapter 10,
• Lecture 16: Driscoll and Braun (2018) Galerkin method,
  
• Today: Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.

17.1 Recap

We are considering linear boundary value problems of the following form:

\begin{align} \mathcal L[u] &:= -u'' + p u' + q u = f &&\text{on } (a,b) \nonumber\\ &u(a)=0, \qquad u(b) = 0 . \tag{BVP} \end{align}

where p, q may be functions of x. By approximating the continuous problem with finite differences, we obtained a discrete problem: find U \in \mathbb R^{n+1} (with index starting at 0) such that

\begin{align} \begin{split} &\mathcal L_h [U] = F \\ &U_0 = U_n = 0 \end{split}\tag{$\text{BVP}_h$} \end{align}

where \mathcal L_h depends on the choice of finite difference scheme. Here is the implementation we saw last time:

"""
    fdm(n, p, q, f; a=0,b=1, A=0,B=0)

Approximate the solution to 
    -u′′ + p u′ + q u = 0 on (a,b)
    u(a) = A, u(b) = B
using central finite differences
"""
function fdm(n, p, q, f; a=0.,b=1., A=0.,B=0.)
    # mesh size
    h = (b-a)/n 
    # interior grid points 
    x = [ a + j *(b-a)/n for j ∈ 1:n-1 ] 

    # Finite Difference Coefficients
    # (using a second order central difference) 
    d = (2/h^2) * ones(n-1) + q.(x)
    d₋  = (-1/h^2)*ones(n-2) - p.(x[2:end])/(2h)
    d₊  = (-1/h^2)*ones(n-2) + p.(x[1:end-1])/(2h)
    Lh = Tridiagonal(d₋, d, d₊)
    F = f.(x)
    F[1] += (1/h^2 + p(x[1])/(2h))*A
    F[end] += (1/h^2 - p(x[end])/(2h))*B
    U = Lh \ F

    # add the boundary conditions
    x = [a, x..., b]
    U = [A, U..., B]

    return x, U
end;

Then, we reformulated the problem as: find u such that

\begin{align} &- (\alpha u')' + \beta u = f \qquad \text{on } (a,b) \nonumber\\ &u(a) = 0, \quad u(b) = 0. \tag{BVP} \end{align}

We assume that \alpha has a continuous derivative, \alpha(x)\geq c>0, \beta is continuous with \beta(x) \geq 0, and f is square integrable: \int_a^b |f|^2 < \infty.

We then saw that we can write down a weak formulation of the problem: find u\in V such that

\begin{align} \mathcal A(u,v) &:= \int_a^b \left[ \alpha(x) u'(x) v'(x) + \beta(x) u(x) v(x) \right] \mathrm{d}x \nonumber\\ % &= \int_a^b f v =: (f,v) \end{align}

for all v \in V. Here, V is a space of functions v satisfying the boundary conditions v(a) = v(b) = 0 and for which we can define \mathcal A(u,v) for all u,v \in V (it takes a bit more work to define exactly what we mean by V - see a course including Sobolev spaces - it is some space that includes smooth funcions satisfying the boundary conditions).

This is called the weak formulation because, if u solves (BVP), then u solves the weak formulation.

The idea behind the Galerkin method is to solve the problem projected to a finite dimensional space: Find u \in V_h such that \mathcal A(u,v) = (f,v) for all v \in V_h := \mathrm{span}\{\phi_1,\dots,\phi_m\}, which leads to the following matrix equation:

\begin{align} A u = F \quad \text{where} \quad &u(x) = \sum_{j=1}^{m} u_j \phi_j(x) \nonumber\\ &A_{ij} := \mathcal A(\phi_i,\phi_j) \nonumber\\ % &F_i := (f, \phi_i) \end{align}

As a result, we have taken the continuous problem (find u\in V such that \mathcal A(u,v) = (f,v) for all v\in V) and approximated it by a discrete problem (find u\in V_h such that \mathcal A(u,v) = (f,v) for all v\in V_h) that can be written as a linear system of equations Au = F. Since we are now in “Linear Algebra World”, we can use results from Chapter 2: Iterative Methods (and last semester Chapter 5: Direct Methods) to solve this system of equations.

Last lecture, we saw an example of (i) a subspace V_h (piecewise linear functions), (ii) building the matrix A (assembly) and (iii) solving the linear system of equations:

"""
    fem(n, α, β, f; a=0, b=1)

Approximate the solution to 
    -(α u′)′ + β u = f on (a,b)
    u(a) = 0, u(b) = 0
using finite element method with hat functions
"""
function fem(n, α, β, f; a=0, b=1)
    h = (b-a) / n
    x = a .+ h .* (0:n)

    α_av = (α.(x[1:n])+α.(x[2:n+1]))/2
    β_av = (β.(x[1:n])+β.(x[2:n+1]))/2
    f_av = (f.(x[1:n])+f.(x[2:n+1]))/2

    K_d  = [ (α_av[i]+α_av[i+1])/h for i ∈ 1:n-1 ]  
    K_dl = [ -α_av[i+1]/h for i ∈ 1:n-2 ]    
    K = SymTridiagonal(K_d, K_dl)

    M_d = [ h/3 *(β_av[i] + β_av[i+1]) for i ∈ 1:n-1 ]
    M_dl = [ h/6 *β_av[i+1] for i ∈ 1:n-2 ]
    M = SymTridiagonal(M_d, M_dl)
    
    F = [ (h/2) * ( f_av[i] + f_av[i+1] ) for i = 1:n-1 ]
   
    u = (K+M) \ F
    u = [0; u; 0]    
    return x, u
end

a, b = 0, 2
α(x) = 1
β(x) = 1
f(x) = -8 + 16x^2 - x^4
u_e(x) = x^2*(4-x^2)

x, u = fem(5, α, β, f; a, b)
x2, u2 = fem(10, α, β, f; a, b)

P = plot(u_e, a, b; 
    title=L"FEM solution to $-u''+u=f$", 
    label="exact solution", linestyle=:dot, color=:black)
plot!(x, u; 
    label=L"n=5", markersize=2, m=2)
plot!(x2, u2; 
    label=L"n=10", markersize=2, m=2)

nn = [2^j for j∈2:12]
err_fem = zeros(length(nn))
for (i,n) ∈ enumerate(nn)
    x, u = fem(n, α, β, f; a, b)
    err_fem[i] = abs.( maximum(u - u_e.(x)) ) 
end

Q = plot( nn, err_fem;
    title="error curves",
    xaxis=:log,yaxis=:log, m=2, xlabel=L"n", ylabel="error",
    label="FEM" )
plot!( nn, nn.^(-2); linestyle=:dot, label=L"O(h^2)")

plot( P, Q, layout=(2,1), size=(500,500) )

Here, we saw that this method converges with rate O(h^2) in the max norm as the step size goes to zero.

Today, we are going to see some of the ideas needed to prove that these methods converge. It will be easier to do this in the so-called energy norm instead of the max-norm:

17.2 Error estimates

This section is mainly based on Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.

Recall, we are approximating the solution to the weak formulation: find u\in V such that

\begin{align} \mathcal A(u,v) = (f,v) \qquad \text{for all } v\in V \tag{weak} \end{align}

with the approximation (projection): find u_h \in V_h such that

\begin{align} \mathcal A(u_h,v_h) = (f,v_h) \qquad \text{for all } v_h\in V_h \tag{G$_h$} \end{align}

(we now put subscript h to distinguish u and u_h and (\text{G}_h) is for “Galerkin” with mesh size h) where

\begin{align} \mathcal A(u,v) &:= \int_a^b \left[ \alpha(x) u'(x) v'(x) + \beta(x) u(x) v(x) \right] \mathrm{d}x \nonumber\\ % (f,v) &:= \int_a^b f v, \end{align}

and \alpha(x) \geq c > 0 and \beta(x) \geq 0.

Under these assumptions on \alpha, \beta, we can define a norm:

Theorem 17.1 The function

\begin{align} \| u \|_{\mathcal A} := \sqrt{\mathcal A(u,u)} \end{align}

defines a norm on V (recall V is a space of functions with zero boundary conditions).

“Proof”.

Positivity: Since \alpha,\beta\geq0, we have

\begin{align} \mathcal A(u,u) % = \int_a^b \left[ \alpha (u')^2 + \beta u^2\right] % \geq 0. \end{align}

Definiteness: Let’s do this for smooth u (as mentioned before, V is a bigger space but I don’t want to introduce Sobolev spaces here). If \mathcal A(u,u) = 0 then \int_a^b \alpha(x) (u'(x))^2\mathrm{d}x = 0, then (since \alpha(x) \geq c>0), we have u' = 0 on [a,b]. Since u\in V, we have u(a) = 0 and so

\begin{align} u(x) = u(a) + \int_a^{x} u'(t) \mathrm{d}t % = 0. \end{align}

That is, \|u\|_{\mathcal A} = 0 implies u = 0.

Homogeneity: Since (\lambda u)'(x) = \lambda u'(x) we have

\begin{align} \| \lambda u \|_{\mathcal A} % &= \sqrt{\mathcal A(u, u)} \nonumber\\ % &= \sqrt{ \int_a^b \left[\alpha(x) (\lambda u'(x))^2 + \beta(x) (\lambda u(x))^2 \right]} \nonumber\\ % &= |\lambda| \|u\|_{\mathcal A}. \end{align}

Triangle inequality: This follows from the fact that \mathcal A is an inner product as we will see: since \mathcal A is symmetric (\mathcal A(u,v) = \mathcal A(v,u)), positive (\mathcal A(u,u) \geq0), and linear in each component (bilinear), we have

\begin{align} 0 \leq \mathcal A(u+\lambda v,u+\lambda v) &= \mathcal A(u,u) + 2\lambda \mathcal A(u,v) + \lambda^2\mathcal A(v,v) \nonumber\\ % &= \|u\|_{\mathcal A}^2 + 2\lambda\mathcal A(u,v) + \lambda^2 \|v\|_{\mathcal A}^2 \end{align}

This is a quadratic in \lambda with at most one real root (since it is positive), and so the discriminant is less than or equal to zero: D := 4\mathcal A(u,v)^2 - 4 \|u\|_{\mathcal A}^2 \|v\|_{\mathcal A}^2 \leq 0. That is, we have the Cauchy-Schwarz inequality \mathcal A(u,v) \leq \|u\|_{\mathcal A} \|v\|_{\mathcal A}. We therefore have (taking \lambda=1)

\begin{align} \mathcal A(u+v,u+v) &= \|u\|_{\mathcal A}^2 + 2\mathcal A(u,v) + \|v\|_{\mathcal A}^2 \nonumber\\ % &\leq \|u\|_{\mathcal A}^2 + 2\|u\|_{\mathcal A}\|v\|_{\mathcal A} + \|v\|_{\mathcal A}^2\nonumber\\ % &= (\|u\|_{\mathcal A} + \|v\|_{\mathcal A})^2, \end{align}

and so we the triangle inequality \|u+v\|_{\mathcal A}\leq \|u\|_{\mathcal A} + \|v\|_{\mathcal A}. \square

It turns out that the solution u_h to (\text{G}_h) is the orthogonal projection of u (the solution to (\text{weak})) to V_h:

Theorem 17.2 (Energy norm) Let u solve \text{(weak)} and u_h solve (\text{G}_h). Then,

(i) Galerkin orthogonality: \mathcal A(u-u_h, v_h) = 0 for all v_h \in V_h.

(ii) \|u-u_h\|_{\mathcal A} = \min_{v_h \in V_h} \|u-v_h\|_{\mathcal A}

Proof.

(i) Using the linearity of \mathcal A and the fact that V_h \subset V, we have

\begin{align} \mathcal A(u - u_h, v_h) &= \mathcal A(u, v_h) - \mathcal A(u_h, v_h)\nonumber\\ % &= (f,v_h) - (f, v_h) = 0 \end{align}

for all v_h \in V_h.

(ii) Using the bilinearity of \mathcal A, we have

\begin{align} \|u-v_h\|_{\mathcal A}^2 &= \mathcal A(u-v_h, u-v_h) \nonumber\\ % &= \mathcal A(u-u_h + u_h - v_h, u-u_h+u_h-v_h) \nonumber\\ % &= \mathcal A(u-u_h, u-u_h) % + 2 \mathcal A(u-u_h, u_h-v_h) \nonumber\\ % &\qquad + \mathcal A(u_h - v_h, u_h-v_h) \nonumber\\ % &= \|u-u_h\|_{\mathcal A}^2 + \|u_h-v_h\|_{\mathcal A}^2. \end{align}

Here, we have used the Galerkin orthogonality \mathcal A(u - u_h, u_h-v_h) = 0. As a result, we have \|u-u_h\|_{\mathcal A} \leq \|u-v_h\|_{\mathcal A} for all v_h \in V_h. \square

This theorem shows that u_h is the best approximation (with respect to the energy norm) to u in the space V_h. We next show how to use this to get an error estimate \| u - u_h \|_{\mathcal A} as a function of the mesh size h.

Recall that the functions \phi_j are cardinal (that is, \phi_j(x_i) = \delta_{ij}). Therefore, the interpolation

\begin{align} I_h u(x) := \sum_{j=1}^m u(x_j) \phi_j(x) \end{align}

is a function that belongs to V_h and satisfies I_h u(x_j) = u(x_j). Therefore, by the theorem above, we have

\begin{align} \| u - u_h \|_{\mathcal A} \leq \| u - I_h u \|_{\mathcal A}. \end{align}

We can obtain error estimates for finite element methods by finding an upper bound on the right hand side of this inequality.

Theorem 17.3 Suppose that u is twice continuously differentiable on [a,b]. Then, piecewise linear finite elements converge with rate O(h) in the energy norm:

\begin{align} \| u - u_h \|_{\mathcal A} \leq C h \| u''\|_{L^2([a,b])} \end{align}

where \| u''\|_{L^2([a,b])} := \sqrt{\int_a^b |u''(t)|^2 \mathrm{d}t} denotes the L^2 norm of u''.

Proof.

We apply Theorem thm-EN and estimate the right hand side of \| u - u_h \|_{\mathcal A}\leq \| u - I_h u \|_{\mathcal A} where I_hu is the interpolation defined above. (The proof that the right hand side of this is O(h) is similar to a proof we had last semester).

We will show that, on defining v := u - I_hu, we have

\begin{align} \|v\|_{L^2([a,b])} &\leq h \|v'\|_{L^2([a,b])} % \leq h^2 \|v''\|_{L^2([a,b])}. \tag{*} \end{align}

If this is true, then since v'' = u'', we have

\begin{align} \| u - u_h \|_{\mathcal A}^2 &\leq \| u - I_hu \|_{\mathcal A}^2 \nonumber\\ % &= \int_a^b \left( \alpha(x) |v'(x)|^2 + \beta(x) |v(x)|^2 \right) \mathrm{d}x \nonumber\\ % &\leq \Big( \max_{[a,b]} |\alpha| \Big) \|v'\|^2_{L^2([a,b])} + \Big( \max_{[a,b]} |\beta| \Big)\|v\|^2_{L^2([a,b])} \nonumber\\ % &\leq \left[ h^2 \Big( \max_{[a,b]} |\alpha| \Big) + h^4 \Big( \max_{[a,b]} |\beta| \Big) \right] \|u''\|^2_{L^2([a,b])}. \end{align}

We conclude by taking square roots.

Therefore, we only need to prove Poincaré’s inequality (*) (Poincaré’s inequality is a bound of the form \|v\|_{L^2} \leq \|v'\|_{L^2} which is true e.g. when v is zero on the boundary, or if v has zero mean \int v = 0). To prove this we need the Cauchy-Schwarz inequaltiy: \int fg \leq \|f\|_{L^2} \|g\|_{L^2} (the proof of this is exactly the same as before:

\begin{align} 0 &\leq (f+\lambda g,f+\lambda g) \nonumber\\ &= \int |f|^2 + 2\lambda \int f g + \lambda^2 \int|g|^2. \end{align}

This is a quadratic in \lambda and the discriminant being negative is equivalent to C-S.)

Recall that v := u - I_hu. Since I_hu is the piecewise linear interpolation of u on \{x_j\}, we have v(x_{k-1}) = v(x_k) = 0. Moreover, since v is twice continuously differentiable, there exists \xi_k\in[x_{k-1},x_k] such that v'(\xi_k) = 0. Therefore, we have

\begin{align} v(x) &= v(x_{k-1}) + \int_{x_{k-1}}^x v'(t) \mathrm{d}t = \int_{x_{k-1}}^x v'(t) \mathrm{d}t \nonumber\\ v'(x) &= v'(\xi_{k}) + \int^{x}_{\xi_k} v''(t) \mathrm{d}t = \int^{x}_{\xi_k} v''(t) \mathrm{d}t. \end{align}

In the following, we will apply the the exact same argument to both w = v and w = v'. On defining y_k := x_{k-1} if w=v and y_k := \xi_{k} if w=v' and apply the Cauchy-Schwarz inequality, we have: for x \in [x_{k-1},x_k],

\begin{align} |w(x)|^2 % &= \left|\int_{y_{k}}^x 1\cdot w'(t) \mathrm{d}t\right|^2 \nonumber\\ % &\leq \left(\int_{y_k}^x 1^2\mathrm{d}t\right) \left(\int_{y_{k}}^x |w'(t)|^2 \mathrm{d}t\right) \nonumber\\ % &\leq h \int_{x_{k-1}}^{x_{k}} |w'(t)|^2 \mathrm{d}t. \end{align}

As a result, we have \int_{x_{k-1}}^{x_k} |w(x)|^2 \mathrm{d}x \leq h^2 \int_{x_{k-1}}^{x_{k}} |w'(t)|^2 \mathrm{d}t and, summing over k, we obtain \|w\|_{L^2([a,b])}^2 \leq h^2 \| w'\|_{L^2([a,b])}^2. \square

17.3 Potential presentation topic

Suli and Mayers (2012) chapter 14.5 - A posteriori error analysis and mesh refinement.

TipTO-DO

• Take a look at the suggestions for Presentations/Posters/Papers and sign up to present,
• The last few lectures have been more difficult: Read through the proofs/exercises from the notes and make sure you are up-to-date,
  
• (Sign-up to office hours if you would like)

Previous reading:

• Chapter 1 reading:
  • Burden, Faires, and Burden (2015) sections 5.1 - 5.6,
    
  • Driscoll and Braun (2018) zero stability
• Chapter 2 reading:
  • Recap matrices: sections 7.1 - 7.2 of Burden, Faires, and Burden (2015)
  • Jacobi & Gauss-Seidel: section 7.3 Burden, Faires, and Burden (2015)
    
  • CG: section 7.6 Burden, Faires, and Burden (2015)
• Chapter 3 reading:
  • Sections 11.1 - 11.2 of Burden, Faires, and Burden (2015): shooting methods
  • Lecture 16: Driscoll and Braun (2018) Galerkin method
  • Today: Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.

Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.

Driscoll, Tobin A, and Richard J Braun. 2018. Fundamentals of Numerical Computation. https://fncbook.com/.

Suli, Endre, and David F Mayers. 2012. An Introduction to Numerical Analysis. Cambridge, England: Cambridge University Press.