16  Galerkin Method and Finite Elements

Lecture 16

Published

April 8, 2026

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

16.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. We will see that including non-zero boundary conditions is easy (as in the case of finite differences). 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.

Today, we will look at an alternative approach involving approximating integrals instead of derivatives.

16.2 Galerkin Method

The content of this lecture is mainly contained in Driscoll and Braun (2018) Galerkin method.

We implement this here (we also have the finite difference code from last lecture):

# this is adapted from [Galerkin method](https://fncbook.com/galerkin/)
"""
    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)

    KK, MM, bb = [1 -1; -1 1], (1/6)*[2 1; 1 2], (1/2)*[1;1]

    α_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, M, b = zeros(n-1,n-1), zeros(n-1,n-1), zeros(n-1)
    K[1,1], K[n-1,n-1] = α_av[1]/h, α_av[n]/h
    M[1,1], M[n-1,n-1] = β_av[1]*h/3, β_av[n]*h/3
    b[1], b[n-1] = f_av[1]*h/2, f_av[n]*h/2

    for k ∈ 2:n-1
        K[k-1:k,k-1:k] += (α_av[k]/h) * KK
        M[k-1:k,k-1:k] += β_av[k]*h * MM
        b[k-1:k] += f_av[k]*h * bb
    end

    u = (K+M) \ b
    u = [0; u; 0]    
    return x, u
end

"""
    fdm(n, p, q, f; a=0,b=1, α=0,β=0)

Approximate the solution to 
    -u′′ + p u′ + q u = 0 on (a,b)
    u(a) = α, u(b) = β
using central finite differences
"""
function fdm(n, p, q, f; a=0.,b=1., α=0.,β=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))*α
    F[end] += (1/h^2 - p(x[end])/(2h))*β
    U = Lh \ F

    # add the boundary conditions
    U = [α, U..., β]
    x = [a, x..., b]

    return x, U
end;

Example 16.1 (Driscoll and Braun (2018) Example 10.6.2) Let’s apply this to the following problem

\begin{align} &-(x^2 u')' + 4u = \sin \pi x && \text{on }(0,1) \nonumber\\ % &u(0) = u(1) = 0 \end{align}

  • Write this in the form -u'' + p u' + q = f,
  • Let’s see how our finite element method compares to the finite difference method we had last time:
p(x) = -2/x
q(x) = 4/x^2
f(x) = sin(π*x)/x^2    

α(x) = x^2
β(x) = 4
g(x) = sin(π*x)

x2, u2 = fem(30, α, β, g)
plot(x2, u2; 
     label="finite element method", markersize=2, m=2)

x, u = fdm(30, p, q, f)
P1 = scatter!(x, u;  
     label="finite difference method", markersize=2)

Here is the error curve for both:

egrid = range(0,1,500)
x, u = fem(2^12, α, β, g)
u_e = linear_interpolation(x, u)

nn = [2^j for j∈2:10]
err = zeros(length(nn))
err2 = zeros(length(nn))

for (i,n) ∈ enumerate(nn)
    x, u = fem(n, α, β, g)
    uf = linear_interpolation(x, u)
    err[i] = abs.( maximum(uf.(egrid) - u_e.(egrid)) ) 

    x, u = fdm(n, p, q, f)
    uf = linear_interpolation(x, u)
    err2[i] = abs.( maximum(uf.(egrid) - u_e.(egrid)) ) 
end

plot( nn, err;
    xaxis=:log,yaxis=:log, m=2, xlabel="n", ylabel="error",
    label="FEM" )
plot!( nn, err;
    label="FDM" )
plot!( nn[3:end], nn[3:end].^(-2), linestyle=:dash, label=L"O(h^2)")

Exercise 16.1 Exercises from Driscoll and Braun (2018)

16.3 Error estimates

Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6 and the idea

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

TipTO-DO
  • Today: Driscoll and Braun (2018) Galerkin method
  • Take a look at the suggestions for Presentations/Posters/Papers and sign up to present,
  • Read through the proofs/exercises from this lecture and Monday,
  • (Sign-up to office hours if you would like)

Previous reading:

  • Chapter 1 reading:
  • 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
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.