13) Introduction to Finite Element discretizations

13) Introduction to Finite Element discretizations#

Last time:

  • Stability of multi-step schemes

Today:

  1. Minimization approach

  2. Variational formulations

  3. Finite Elements

  4. Quadrature

  5. Gauss-Legendre quadrature

  6. Lagrange Interpolating Polynomials

  7. Galerkin Spectral Element Method (SEM): \(L^2\) projection

1. Minimization approach#

Let’s consider the Poisson’s problem once again

\[\begin{split} \begin{align*} - \Delta u &= f \qquad \textrm{ in } \Omega \qquad {\tag{13.1}}\\ u &= g \qquad \textrm{ on } \partial \Omega \end{align*} \end{split}\]

This is for the diffusivity constant \(\kappa = 1\), but we will use its most general form later.

Recall that the general theory of solutions to Laplace’s equation is known as potential theory. In particular, using what is called Riesz potential.

In physics, the solution to Poisson’s equation is the potential field caused by a given electric charge or mass density distribution; with the potential field known, one can then calculate the corresponding electrostatic or gravitational (force) field. One way to solve this is by the minimization approach, justified by the fact that the solution will exhibit minimal energy (since it’s a steady state problem).

Define the total energy as

\[ \varepsilon(u) = \int_{\Omega} \left( \frac{1}{2} |\nabla u|^2 -fu \right) dx \]

Then

\[ u = \arg\min \varepsilon (w ) \quad, \, w \in C^2(\Omega) \cap C(\bar{\Omega}) , \, w = g \textrm{ on } \partial \Omega \]

How do we show that \(u\) satisfies the minimal energy requirement? Let’s perturb the equilibrium of a small quantity and let’s observe the behavior.

Let \(w = u + v\), (consider \(v = du\), a small perturbation of \(u\)) with \(v = 0\) on \(\partial \Omega\). Then

\[\begin{split} \begin{align*} \varepsilon(w) &= \int_{\Omega} \left( \frac{1}{2} |\nabla u^2 + \nabla v^2| -fu -fv \right) dx \\ &= \varepsilon(u) + \frac{1}{2}\int_{\Omega} | \nabla v^2| x + \int_{\Omega} \left( \nabla u \nabla v -fv \right) dx \end{align*} \end{split}\]

Consider

\[ 0 = \int_\Omega \nabla \cdot (v \nabla u)= \int_{\Omega} \left( \nabla u \nabla v + \underbrace{\cancel{ \Delta u v}}_{= -fv\textrm{, b/c u is sol}} \right) dx = \int_{\Omega} \left( \nabla u \nabla v -fv \right) dx \qquad \tag{13.2} \]

Hence the last integral on the right cancels and we have that the minimal energy was an equilibrium solution.

  • Considering a small perturbation \(v = \delta u\), so that \(w = u + \delta u\), is the reason why much of this theory is based on what is called the variational formulation (as \( \delta v\) is a small variation from \(u\)).

The right-most integral

\[ \int_{\Omega} \left( \nabla u \nabla v -fv \right) dx = 0 \]

gives

\[ \int_{\Omega} \nabla u \nabla v dx = \int_{\Omega} fv dx \tag{13.3} \]

which is the variational formulation or weak form of (13.1) that, in contrast, is called the strong form.

2. Variational formulations#

We have just derived this variational formulation from an energy minimization argument. But, more simply, we can think of obtaining it equivalently by considering the integral form of our PDE and applying integration by parts and its variations in multiple dimensions, which relate to Gauss’ divergence theorem.

From Calculus we know that:

If \(f(x)\) is a minimum, then \(df = f'(x) dx = 0\) for all \(dx\).

So if \(u \in \mathcal V = \{ w \in C^2(\bar{\Omega}) | w = 0 \textrm { on } \partial \Omega \}\) is a minimizer, then

\[ \begin{align*} \int_\Omega \underbrace{\lVert \nabla u \rVert^{p-2}}_{\kappa(u)} \nabla u \cdot \nabla du = 0 \end{align*} \]

for all \(du \in \mathcal V\).

  • It’s common to use \(v\) in place of \(du\).

Integration by parts#

One dimension#

it’s the product rule backwards

\[\begin{split} \begin{align*} \int_a^b d(uv) &= \int_a^b u dv + \int_a^b v du \\ (uv)_a^b &= \int_a^b u dv + \int_a^b v du \\ \int_a^b u dv &= (uv)_a^b - \int_a^b v du \end{align*} \end{split}\]

  • You can move the derivative to the other term; it’ll cost you a minus sign and a boundary term

Multiple dimensions#

\[ \begin{align*} \int_\Omega v \nabla\cdot \mathbf f = -\int_\Omega \nabla v \cdot \mathbf f + \int_{\partial \Omega} v \mathbf f \cdot \mathbf n \end{align*} \]

Definition: Weak Derivative

Let \(u\) be a function in the Lebesgue space \(L^1([a,b])\). We say that \(v\) in \(L^1([a,b])\) is a weak derivative of \(u\) if

\[ \int_a^b u(x) \phi '(x) dx = - \int_a^b v(x) \phi(x) dx \]

for all inifinitely differentiable functions \(\phi\) such that \(\phi(a)=\phi(b)=0\).

Generalizing in \(n\) dimensions, if \(u\) and \(v\) are in the space of locally integrable functions \(L^1_{\textrm{Loc}}(\Omega)\) for some open set \(\Omega \subset \mathbb{R}^n\), and if \(\alpha\) is a multi-index, we say that \(v\) is the \(\alpha\)-th weak derivative of \(u\) if

\[ \int_{\Omega} u D^{\alpha} \phi = (-1)^{|\alpha|} \int_\Omega v \phi \]

for all \(\phi \in C^{\infty}(\Omega)\), that is, for all inifinitely differentiable functions \(\phi\) with compact support in \(\Omega\). Here \(D^{\alpha}\phi\) is defined as

\[ D^\alpha \phi = \frac{\partial^{|\alpha|} \phi}{\partial x_1^{\alpha_1} \cdots \partial x_n^{\alpha_n}} \]

The basic recipe for turning a PDE into a variational problem is:

  • Multiply the PDE by a function \(v\)

  • Integrate the resulting equation over the domain \(\Omega\)

  • Perform integration by parts of those terms with second order derivatives

The function \(v\) which multiplies the PDE is called a test function. The unknown function \(u\) that is to be approximated is referred to as a trial function.

We have already seen the variational form or weak form of the Poisson’s equation (13.3) above,

\[ \int_{\Omega} \nabla u \nabla v dx = \int_{\Omega} fv dx \tag{13.3} \]

  • A rule of thumb when deriving variational formulations is that one tries to keep the order of derivatives of \(u\) and \(v\) as small as possible.

  • Another feature of variational formulations is that the test function \(v\) is required to vanish on the parts of the boundary where the solution \(u\) is known. This means that \(v\) is \(0\) on the whole boundary \(\partial\Omega\).

If we require that equation (13.3) holds for all test functions \(v\) in some suitable space \(\hat{\mathcal{V}}\), the so-called test space, we obtain a well-defined mathematical problem that uniquely determines the solution \(u\) which lies in some function space \(\mathcal{V}\). Note that \(\mathcal{V}\) does not have to be the same space as \(\hat{\mathcal {V}}\). We call the space \(\mathcal{V}\) the trial space.

The proper way to state the variational problem is:

  • Find \(u\in \mathcal{V}\) such that

\[\int_\Omega \nabla u \cdot \nabla v~\mathrm{d} x = \int_\Omega f v~\mathrm{d} x\qquad \forall v \in \hat{\mathcal{V}}.\]

For the present problem, the trial and test spaces \(\mathcal{V}\) and \(\hat{\mathcal{V}}\) are defined as

\[\begin{split}\begin{align*} \mathcal{V} &=\{v\in H^1(\Omega) \vert v=u_D&&\text{on } \partial \Omega \},\\ \hat{\mathcal{V}}&=\{v\in H^1(\Omega) \vert v=0 &&\text{on } \partial \Omega \}. \end{align*} \end{split}\]

In short, \(H^1(\Omega)\) is the Sobolev space containing functions \(v\) such that \(v^2\) and \(\vert \nabla v \vert ^2\) have finite integrals over \(\Omega\).

The solution of the underlying PDE must lie in a function space where the derivatives are also continuous, but the Sobolev space \(H^1(\Omega)\) allows functions with discontinuous derivatives.

This weaker continuity requirement in our weak formulation (caused by the integration by parts) is of great importance when it comes to constructing the finite element function space. In particular, it allows the use of piecewise polynomial function spaces. This means that the function spaces are constructed by stitching together polynomial functions on simple subdomains such as intervals, triangles, quadrilaterals, tetrahedra and hexahedra.

The variational problem is a continuous problem: it defines the solution \(u\) in the infinite-dimensional function space \(\mathcal{V}\).

The finite element method for the Poisson equation finds an approximate solution of the variational problem by replacing the infinite-dimensional function spaces \(\mathcal{V}\) and \(\hat{\mathcal{V}}\) by discrete (finite dimensional) trial and test spaces \(\mathcal{V}_h\subset \mathcal{V}\) and \(\hat{\mathcal{V}}_h \subset \hat{\mathcal{V}}\).

The discrete variational problem reads: Find \(u_h\in \mathcal{V}_h\) such that

\[ \begin{align*} \int_\Omega \nabla u_h \cdot \nabla v~\mathrm{d} x &= \int_\Omega fv~\mathrm{d} x && \forall v \in \hat{\mathcal{V}}_h. \end{align*} \]

This variational problem, together with suitable definitions of \(\mathcal{V}_h\) and \(\hat{\mathcal{V}}_h\) uniquely define our approximate numerical solution of the Poisson equation.

Note that the boundary condition is encoded as part of the test and trial spaces. This might seem complicated at first glance, but means that the finite element variational problem and the continuous variational problem look the same.

3. Finite Elements#

Definition: Finite Element (Ciarlet 1978 definition)

A finite element is a triple \((K, P, N)\) where

  • \(K\) is a bounded domain with piecewise smooth boundary

  • \(P = \{p_j\}_1^k\) is a finite-dimensional function space over \(K\) (“prime basis”)

  • \(N = \{n_i\}_1^k\) is a basis for the dual space (“nodes”)

It is common for nodes to be pointwise evaluation

\[ \langle n_i, p_j \rangle = p_j(x_i) \]

We will introduce the following notation for variational problems:

  • Find \(u\in \mathcal{V}\) such that

\[ \begin{align*} a(u,v)&=L(v)&& \forall v \in \hat{\mathcal{V}}. \end{align*} \]

For the Poisson equation, we have:

\[\begin{split} \begin{align*} a(u,v) &= \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d} x,\\ L(v) &= \int_{\Omega} fv~\mathrm{d} x. \end{align*} \end{split}\]

In the literature \(a(u,v)\) is known as the bilinear form and \(L(v)\) as a linear form.

  • In general, for every linear problem, we want to identify all terms with the unknown \(u\) and collect them in \(a(u,v)\), and collect all terms with only known functions in \(L(v)\).

We note that the bilinear form \(a(u,v)\) is

  • bilinear (linear in each of its arguments)

  • symmetric: \(a(u, v) = a(v,u)\)

  • positive definite: \(a(u, u) > 0\) when \(u \ne 0\) thus defines an inner product on the function space \(\mathcal{V}\).

We also introduce the \(L^2\) inner product

\[ \langle u, v \rangle = \int_\Omega u(x) v(x) ~\mathrm{d} x \]
so that our continuous weak form variational problem is to find \(u \in \mathcal{V}\) such that
\[ a(v, u) = \langle v, f \rangle, \quad \forall v\in \mathcal{V}. \]

The discrete variational problem, also known as Galerkin discretization, is to find \(u_h \in V_h \subset \mathcal{V}\) such that

\[ a(v_h, u_h) = \langle v_h, f \rangle, \quad \forall v_h \in \mathcal{V}_h . \]
Since \(\mathcal{V}_h \subset \mathcal{V}\), we can subtract these two, yielding
\[ a(v_h, u_h - u) = 0, \quad \forall v_h \in \mathcal{V}_h .\]
This says that the error in the discrete solution \(u_h - u\) is \(a\)-orthogonal to all test functions \(v_h\).

Galerkin optimality via energy norms#

We can also define the “energy norm” or \(a\)-norm,

\[ \lVert u \rVert_a = \sqrt{a(u,u)} . \]
This norm satisfies the Cauchy-Schwarz inequality,
\[ \lvert a(u,v) \rvert \le \lVert u \rVert_a \lVert v \rVert_a . \]
Now,

\[\begin{split} \begin{align*} \lVert u_h - u \rVert_a^2 &= a(u_h - u, u_h - u) \\ &= a(u_h - v_h, u_h - u) + a(v_h - u, u_h - u) \\ &= a(v_h - u, u_h - u) \\ &\le \lVert v_h - u \rVert_a \lVert u_h - u \rVert_a . \end{align*} \end{split}\]

In other words,

\[ \lVert u_h - u \rVert_a \le \lVert v_h - u \rVert_a, \quad \forall v_h \in \mathcal{V}_h . \]

So the solution \(u_h\) computed by the Galerkin discretization is optimal over the subspace \(\mathcal{V}_h\) as measured in the \(a\)-norm.

Observations#

  • The Galerkin method computes the exact solution any time it resides in the subspace \(\mathcal{V}_h\).

  • The Galerkin method is automatically symmetric any time the weak form is symmetric.

  • The Galerkin method can be spectrally (i.e., exponentially) accurate, similar to the Chebyshev finite difference methods.

  • For a nonlinear problem, discretization and differentiation will commute.

4. Quadrature#

We have seen that the variational form is an integral form of the PDE. Therefore, in Finite Element Methods (FEM), we need a way to compute integrals. We will do this via numerical quadrature rules.

For a \(1\)-D domain, our elements will be a line segment and we’ll define quadrature on the standard line segment, \(\hat \Omega = (-1, 1)\).

A quadrature is a set of points \(q_i\) and weights \(w_i\) such that

\[ \int_{-1}^1 f(x) dx \approx \sum_{i=1}^n w_i f(q_i) = \mathbf{w}^T f(\mathbf{x}) \]

Integration via polynomial interpolation#

  • Pick some points \(x_i\)

  • Evaluate the function there \(f_i = f(x_i)\)

  • Find the polynomial that interpolates the data

  • Integrate the polynomial

Questions:

  • What order of accuracy can we expect?

  • What degree polynomials can be integrate exactly?

Answer:

  • \(n\)-point Gauss quadrature exactly integrates polynomials of degree \(2n-1\)

Legendre polynomials#

The Legendre polynomials \(P_0(x) = 1\), \(P_1(x) = x\), …, are pairwise orthogonal

\[\int_{-1}^1 P_m(x) P_n(x) = 0, \quad \forall m\ne n.\]

They can be used as a set of basis for Gauss quadrature, giving rise to Gauss-Legendre quadrature.

using LinearAlgebra
using Plots
using LaTeXStrings
default(linewidth=4, legendfontsize=12)

# We can define the interpolation problem using the Legendre polynomials as basis

function vander_legendre(x, k=nothing)
    if isnothing(k)
        k = length(x) # Square by default
    end
    m = length(x)
    Q = ones(m, k)
    Q[:, 2] = x
    for n in 1:k-2
        Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1) # definition of Legendre's polynomials via the generating function
    end
    Q
end

vander_legendre (generic function with 2 methods)

x = LinRange(-1, 1, 100)
P = vander_legendre(x, 10)
plot(x, P[:,end], legend=:none) # the last of ten Legendre polynomials

x = LinRange(-1, 1, 50)
P = vander_legendre(x, 6)
plot(x, P) # the first six Legendre polynomials

# Chebyshev nodes
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

function quad_legendre(a, b; n=20)
    x = CosRange(-1, 1, n)
    P = vander_legendre(x)
    x_ab = (a+b)/2 .+ (b-a)/2*x # change of variables for integrals [a,b] -> [-1,1]
    w = (b - a) * inv(P)[1,:]
    x_ab, w
end

function fint_legendre(f, a, b; n=20)
    x, w = quad_legendre(a, b, n=n)
    w' * f.(x)
end

fint_legendre(x -> 1 + x, -1, 1, n=4)

2.0

x, w = quad_legendre(-1, 1)

w' * cos.(x) - fint_legendre(cos, -1, 1)

0.0

How accurate is it?

function plot_accuracy(fint, tests, ns; ref=[1,2])
    a, b = -2, 2
    p = plot(xscale=:log10, yscale=:log10, xlabel="h", ylabel="error",
        legend=:bottomright)
    hs = (b - a) ./ ns
    for (f, F) in tests
        Is = [fint(f, a, b, n=n) for n in ns]
        Errors = abs.(Is .- (F(b) - F(a)))
        scatter!(hs, Errors, label=f)
    end
    for k in ref
        plot!(hs, hs.^k, label="\$h^{$k}\$")
    end
    p
end

plot_accuracy (generic function with 1 method)

# a couple of tests
F_expx(x) = exp(2x) / (1 + x^2)
f_expx(x) = 2*exp(2x) / (1 + x^2) - 2x*exp(2x)/(1 + x^2)^2

F_dtanh(x) = tanh(x)
f_dtanh(x) = cosh(x)^-2

integrands = [f_expx, f_dtanh]
antiderivatives = [F_expx, F_dtanh]
tests = zip(integrands, antiderivatives)

zip(Function[f_expx, f_dtanh], Function[F_expx, F_dtanh])

p = plot_accuracy(fint_legendre, tests, 4:50, ref=1:5) # in a log-log scale

  • It’s spectral (or exponential) convergence!

  • This spectral convergence (which is given by using this specific set of Legendre polynomials as basis and the nonuniform Chebyshev nodes as quadrature points) is the reason why we call this method the Spectral Element Method (SEM), which is a special case of Finite Element Method (FEM).

Prime to nodal basis#

A nodal basis \(\{ \phi_j(x) \}\) is one that satisfies \( n_i(\phi_j) = \delta_{ij} . \)

We write \(\phi_j\) in the prime basis by solving with the generalized Vandermonde matrix \(V_{ij} = \langle n_i, p_j \rangle\),

\[ \phi_j(x) = \sum_k p_k(x) (V^{-1})_{k,j} . \]

k = 10
xn = LinRange(-1, 1, k)
V = vander_legendre(xn)
xx = LinRange(-1, 1, 50)
Pxx = vander_legendre(xx, k)
Pn = Pxx / V # find the basis ϕ as p(x) * inv(V)
plot(xx, Pn)
scatter!(xn, one.(xn))

5. Gauss-Legendre quadrature#

If you were to do it step-by-step, it is a pretty involved numerical process:

  1. Solve for the points, compute the weights

  • Use a Newton solver to find the roots. You can use the recurrence to write a recurrence for the derivatives.

  • Create a Vandermonde matrix and extract the first row of the inverse or (using more structure) the derivatives at the quadrature points.

  1. Use duality of polynomial roots and matrix eigenvalues.

Luckily, packages/libraries come to the rescue.

We want to approximate

\[ \int_{-1}^{-1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i) \]

Note:

To approximate integrals on a general interval \([a,b]\), the problem needs to be translated to \([-1,1]\) with a simple change of variables. Using the substitution \(t = (2x - a - b) / (b-a)\), we find

\[ \int_a^b f(x) dx = \int_{-1}^{1} f \left( \frac{(b-a)t + b + a}{2} \right) \frac{(b-a)}{2} dt \]

In the above, \(w_i\) are the weights defined as

\[ w_i = \int_{-1}^{1}L_i(x) dx \quad i=1, \dots, n \]

with \(L_i\) the Lagrange interpolating polynomial (whose definition is given below) with \(n\) interpolating nodes, and the quadrature nodes \(x_i\) are the roots of the \(n\)-th Legendre polynomial on \([-1,1]\).

function fint_gauss(f, a, b; n=4)
    """Gauss-Legendre integration using Golub-Welsch algorithm"""
    beta = @. .5 / sqrt(1 - (2 * (1:n-1))^(-2))
    T = diagm(-1 => beta, 1 => beta)
    D, V = eigen(T)
    w = V[1,:].^2 * (b-a)
    x = (a+b)/2 .+ (b-a)/2 * D
    w' * f.(x)
end
fint_gauss(sin, -2, 3, n=4)

0.5733948071694299

plot_accuracy(fint_gauss, tests, 3:30, ref=1:4) # in a log-log scale

6. Lagrange Interpolating Polynomials#

Suppose we are given function values \(y_0, \dotsc, y_m\) at the distinct points \(x_0, \dotsc, x_m\) and we would like to build a polynomial of degree \(m\) that goes through all these points. This explicit construction is attributed to Lagrange (though it is said that he was not first):

\[ l_i(x)= \frac{(x-x_0)(x-x_1) \cdots (x-x_{i-1})(x-x_{i+1}) \cdots (x-x_m)}{(x_i-x_0)(x_i-x_1) \cdots (x_i-x_{i-1})(x_i-x_{i+1}) \cdots (x_i-x_m)} \qquad , \; i=0, 1, \ldots, \, m \ . \]

By construction, this is a degree \(m\) polynomial, with \(l_i(x_j) = \delta_{ij}\) (where \(\delta_{ij}\) denotes the Kronecker delta). This just says that \(l_i(x)\) has zeros at all the \(x_j\)’s except \(x_i\), where it is identically equal to \(1\). Hence, we notice that the Lagrange basis for polynomials of degree \(≤ m\) for those nodes is the set of polynomials \(\{l_0(x), l_1(x), \ldots, l_m(x) \}\) which take values \(l_j(x_k) = 0\) if \(k \neq j\) and \(l_j(x_j)=1\).

Lagrange’s interpolating formula

\[\begin{split} p(x) = \sum_{i=0}^m y_i l_i(x) \equiv \sum_{i=0}^m y_i \prod^{m}_{\substack{i = 0\\i \ne j}} \frac{x - x_j}{x_i - x_j} \end{split}\]

FastGaussQuadrature.jl#

Instead of computing quadrature weights and nodes by hand, some useful packages can help us.

using Pkg
Pkg.add("FastGaussQuadrature")
using FastGaussQuadrature

n = 14
x, q = gausslegendre(n)
scatter(x, q, label="Gauss-Legendre", ylabel="weight", xlims=(-1.1, 1.1))
scatter!(gausslobatto(n)..., label="Gauss-Lobatto")

    Updating registry at `~/.julia/registries/General.toml`

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

┌ Warning: Circular dependency detected.
│ Precompilation will be skipped for dependencies in this cycle:
│  ┌ Symbolics → SymbolicsForwardDiffExt
│  └─ Symbolics → SymbolicsPreallocationToolsExt
└ @ Pkg.API.Precompilation ~/julia-1.10.10/share/julia/stdlib/v1.10/Pkg/src/precompilation.jl:583

k = 10
xn = LinRange(-1, 1, k) # using equispaced points
V = vander_legendre(xn) # gauss-Legendre on equispaced points
Pn = Pxx / V
plot(xx, Pn)
scatter!(x, one.(xn))

xn, _ = gausslobatto(k)
V = vander_legendre(xn) # Gauss-Legendre on Lobatto points
Pn = Pxx / V
plot(xx, Pn)
scatter!(xn, one.(xn))

  • Gauss-Lobatto points are preferred in Continuous Galerkin (CG) formulations, that is, where continuity of the solution is desired.

  • On the other hand, for problems in which we have discontinuity of the solution (and we want to be able to capture that, for instance in shock-capturing problems), as in the Discontinuous Galerkin (DG) formulations, Gauss-Lengendre points can be used.

function vander_legendre_deriv(x, k=nothing)
    if isnothing(k)
        k = length(x) # Square by default
    end
    m = length(x)
    Q = ones(m, k)
    dQ = zeros(m, k)
    Q[:, 2] = x
    dQ[:, 2] .= 1
    for n in 1:k-2
        Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
        dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]
    end
    Q, dQ
end

vander_legendre_deriv (generic function with 2 methods)

6. Galerkin Spectral Element Method (SEM): \(L^2\) projection#

A nice test problem that doesn’t require derivatives or boundary conditions: Find \(u \in \mathcal V_h\) such that

\[ \int_{-1}^1 v(x) \big[ u(x) - f(x) \big] = 0, \quad \forall v \in \mathcal{V}_h\]

This is typically called the Mass matrix or Mass operator problem, because of the matrix that arises on the left-hand side. Let’s see its details here.

The Mass operator is defined via the \(L^2\) projection problem, posed as a weak form on a Hilbert space \(\mathcal{V}^p \subset H^1\), i.e., find \(u \in \mathcal{V}^p\) such that for all \(v \in \mathcal{V}^p\)

\[ \langle v,u \rangle = \langle v,f \rangle , \]

where \(\langle v,u\rangle\) and \(\langle v,f\rangle\) express the continuous linear forms defined on \(\mathcal{V}^p\), and, for sufficiently regular \(u\), \(v\), and \(f\), we have:

\[\begin{split} \begin{aligned} \langle v,u \rangle &:= \int_{\Omega} \, v \, u \, dx ,\\ \langle v,f \rangle &:= \int_{\Omega} \, v \, f \, dx . \end{aligned} \end{split}\]

Following the standard finite/spectral element approach, we formally expand all functions in terms of basis functions, such as

\[\begin{split} \begin{align} u(\mathbf x) &= \sum_{j=1}^p u_j \, \phi_j(\mathbf x) ,\\ v(\mathbf x) &= \sum_{i=1}^p v_i \, \phi_i(\mathbf x) . \end{align} \end{split}\]

The coefficients \(\{u_j\}\) and \(\{v_i\}\) are the nodal values of \(u\) and \(v\), respectively, and we see that \(\mathcal{V}^p = \textrm{span}\{ \phi_1, \phi_2,\ldots , \phi_p \}\). Inserting these expressions in the weak form, we obtain the inner products:

\[ \langle v,u \rangle = \mathbf v^T M \mathbf u , \qquad \langle v,f\rangle = \mathbf v^T \mathbf b \,. \]

Here, we have introduced the mass matrix, \(M\), and the right-hand side, \(\mathbf b\),

\[ M_{ij} := (\phi_i,\phi_j), \;\; \qquad b_{i} := \langle \phi_i, f \rangle, \]

each defined for index sets \(i,j \; \in \; \{1,\dots,p\}\), where \(p\) is the degree of the polynomial basis chosen (defined on \(P = p+1\) points).

function febasis(P, Q)
    x, _ = gausslobatto(P)
    q, w = gausslegendre(Q)
    Pk, _ = vander_legendre_deriv(x)
    Bp, Dp = vander_legendre_deriv(q, P)
    B = Bp / Pk
    D = Dp / Pk
    x, q, w, B, D
end

function L2_galerkin(P, Q, f)
    x, q, w, B, _ = febasis(P, Q)
    A = B' * diagm(w) * B
    rhs = B' * diagm(w) * f.(q)
    u = A \ rhs
    x, u
end

L2_galerkin (generic function with 1 method)

f(x) = tanh(5*x)
P, Q = 10, 30
x, u = L2_galerkin(P, Q, f)
error = u - f.(x)
plot(x, u, marker=:auto, legend=:bottomright, label = "approximate")
plot!(f, title="Galerkin error: $(norm(error))", label = L"\tanh(5x)")

Observations

  • Being the \(\{ \phi_1, \phi_2,\ldots , \phi_p \}\) basis functions a basis for the space \(\mathcal{V}^p\), they are linearly independent, hence each of them is orthogonal to the others.

  • This makes the Mass matrix \(M\) full-rank and invertible.

Embed the solution in the continuous space#

In thel last figure, the Galerkin solution was represented by a piecewise linear function, emphasizing its nodal values. Let’s make a nicer (smoother) plot by multiplying by the corresponding nodal basis function.

P, Q = 18, 18
x, u = L2_galerkin(P, Q, f)
error = u - f.(x)
_, xx, _, BB, _ = febasis(P, 100)
plot(f, label="exact", title="Galerkin error: $(norm(error))",
  xlim=(-1, 1), legend=:bottomright)
plot!(xx, BB * u, label="Galerkin solution")
scatter!(x, u, marker=:circle, color=:black, label="Nodal values")

How accurate is it?#

ns = 2 .^ (2:8)
function L2_error(n)
    x, u = L2_galerkin(n, n, f)
    norm(u - f.(x))
end
plot(ns, L2_error.(ns), marker=:circle, yscale=:log10, xscale=:log10,
  title="Convergence", xlabel="number of points", ylabel="error", label = "L2 Galerkin error")
plot!([n -> n^-2, n -> n^-4], label=["\$n^{-2}\$" "\$n^{-4}\$"])
Contents