Galerkin methods

In the previous section, we observed that we can approximate a known function \(g\) by a set of piecewise polynomial functions \(\phi_i\).

We did not go much detail about how these functions would look like. In this section, we look at this is some more detail.

We start by looking at an interval, subdivided into a set of elements. For each of these segments, we will define a local basis function and a set of global basis functions.

Composition of \(u_h\)

Given a finite dimensional function space \(V_h\), with basis functions \(\phi_i, i=0,\dots,N-1\), we can write any function \(u_h\) in \(V_h\) as

\[ \begin{align*} u_h(x) &= \sum_{i=0}^{N-1} u_i \phi_i(x) \end{align*} \]

where \(u_i\) are the coefficients of function. These coefficients are also known as the degrees of freedom DOFs of the function.

Choice of basis functions

We can do the same thing for a higher order set of polynomials

However, in the setting of a PDE, we do not know the solution \(u\).

In this section, we will explain how to find the solution \(u\) to a PDE using the Galerkin method.

We will consider the Poisson equation in 1D as a starting point

\[\begin{split} \begin{align*} -\frac{\partial^2 u}{\partial x^2} &= f(x) \quad \text{in} \quad \Omega = [0, 1] \\ u(0) &=0 \\ u(1) &=0 \end{align*} \end{split}\]

We rewrite the equation above as a residual

\[ \begin{align*} R(u) = -\frac{\partial^2 u}{\partial x^2} - f \end{align*} \]

Goal is to find a \(u\) such that \(R(u)=0\).

If we substitute \(u_h\) into the Poisson equation, we get

\[ \begin{align*} -\sum_{i=0}^{N-1}u_i\frac{\partial^2 \phi_i}{\partial x^2} - f(x) &= 0 \quad \text{in} \quad \Omega = [0, 1] \end{align*} \]

Does this equation make sense for piecewise linear polynomials?

As the double derivative of a linear function is \(0\), we cannot directly solve this equation!

We want to weaken the requirement of the equality

Inner product:

\[ \begin{align*} (a, b) &= \int_\Omega a \cdot b ~\mathrm{d}x \end{align*} \]

With this, we want to find a \(u_h\) such that

\[ \begin{align*} (R(u_h), v) = (0, v) &= 0 \qquad \forall v \in V_h \end{align*} \]

Reconsidering the equation above, we have

\[ \begin{align*} \int_\Omega -\frac{\partial^2 u_h}{\partial x^2}\cdot v ~\mathrm{d} x = \int_{\Omega} f\cdot v~\mathrm{d}x \end{align*} \]

Integrate by parts

\[ \begin{align*} \int_\Omega \frac{\partial u_h}{\partial x}\cdot \frac{\partial v}{\partial x}~\mathrm{d} x - \left[\frac{\partial u_h}{\partial x}v\right]_{0}^{1} = \int_{\Omega} f\cdot v~\mathrm{d}x \end{align*} \]

The bilinear and linear form

We define the bilinear form \(a(u,v)\)

\[ a(u,v)= \int_\Omega \frac{\partial u_h}{\partial x}\cdot \frac{\partial v}{\partial x}~\mathrm{d} x - \left[\frac{\partial u_h}{\partial x}v\right]_{0}^{1} \]

and the linear form \(L(v)\)

\[ L(v) = \int_{\Omega} f\cdot v~\mathrm{d}x \]

We can now insert our expansion of \(u_h\) (where each \(u_i\) is to be determined)

\[ \begin{align*} \sum_{i=0}^{N-1}u_i \int_\Omega \frac{\partial \phi_i}{\partial x}\frac{\partial v}{\partial x}~\mathrm{d} x - \sum_{i=0}^{N-1}u_i\left[\frac{\partial \phi_i}{\partial x}v]\right]_{0}^{1} = \int_{\Omega} f(x)v~\mathrm{d}x \end{align*} \]

By choosing \(v=\phi_0, \dots, \phi_{N-1}\), we can set up a linear system

\[ \begin{align*} Au = b \end{align*} \]

where

\[ \begin{align*} A_{j, i} = \int_\Omega \frac{\partial \phi_i}{\partial x}\frac{\partial \phi_j}{\partial x}~\mathrm{d} x - \sum_{i=0}^{N-1}u_i\left[\frac{\partial \phi_i}{\partial x}\phi_j]\right]_{0}^{1} \end{align*} \]

\[ \begin{align*} b_{j} = \int_{\Omega} f \phi_j~\mathrm{d}x \end{align*} \]

Galerkin orthogonality

A key property of this method is that if \(a(u_h, v) = L(v)\), then the error between the exact solution \(u\) and \(u_h\), given as \(\epsilon=u-u_h\) has the following property:

\[ a(\epsilon, v) = a(u-u_h, v) = a(u, v) - a(u_h, v) = L(v) - L(v) = 0 \]

Thus the error is orthogonal to the chosen sub-space

Questions

How does the matrix \(A\) look like? Is it dense or sparse?

The matrix will be sparse, as the support (function being non-zero) of \(\phi_i\) and \(\phi_j\) is localized to adjacent elements.

Does the Galerkin orthogonality property say something about the quality of the solution?

The Galerkin orthogonality states that you get the best approximation in the chosen function space. However, with a poor choice of function space, the error can be large.

How can we compute the integral \(\int_\Omega \frac{\partial \phi_i}{\partial x}\frac{\partial \phi_j}{\partial x}~\mathrm{d} x\)?

We split the domain \(\Omega\) into non-overlapping elements \(k\in \mathcal{K}\) such that \(\cup_{k\in\mathcal{K}} k = \Omega\) where \(\phi_i\) and \(\phi_j\) are continuous on each section, and use that \(\int_{\Omega} h~\mathrm{d}x = \sum_{k\in\mathcal{K}} \int_{k}h~\mathrm{d}x\).