The Galerkin method - Fundamentals of Numerical Computation

Using finite differences we defined a collocation method in which an approximation of the differential equation is required to hold at a finite set of nodes. In this section we present an alternative based on integration rather than differentiation. Our presentation will be limited to the linear BVP

u'' = p(x)\,u' + q(x)\,u + r(x), \quad a \le x \le b, \quad u(a)=0,\; u(b)=0.

(10.6.1)

However, we will assume that the linear problem is presented in the equivalent form

- \frac{d }{d x} \Bigl[ c(x)\, u'(x) \Bigr] + s(x) \, u(x) = f(x), \quad u(a)=0,\; u(b)=0.

(10.6.2)

Such a transformation is always possible, at least in principle (see Exercise 10.6.5). As with finite differences, a nonlinear problem is typically solved by using a Newton iteration to create a sequence of linear problems.

10.6.1Weak formulation¶

Let (10.6.2) be multiplied by a generic function $\psi(x)$ called the test function, and then integrate both sides in $x$ :

\begin{split} \int_a^b f(x)\psi(x) \,dx &= \int_a^b \bigl[ -(c(x)u'(x))'\psi(x) + s(x)u(x)\psi(x) \bigr] \,dx \\ &= \Bigl[-c(x)u'(x)\psi(x) \Bigr]_{\,a}^{\,b} + \int_a^b \bigl[ c(x)u'(x)\psi'(x) + s(x)u(x)\psi(x)\bigr] \, dx. \end{split}

(10.6.3)

The last line above used an integration by parts.

We now make an important and convenient assumption about the test function. The first term in (10.6.3), consisting of boundary evaluations, disappears if we require that $\psi(a)=\psi(b)=0$ . Doing so leads to

\int_a^b \bigl[ c(x)u'(x)\psi'(x) + s(x)u(x)\psi(x)\bigr] \,dx = \int_a^b f(x)\psi(x) \,dx,

(10.6.4)

which is known as the weak form of the differential equation (10.6.2).

Every solution of (10.6.2) (what we might now call the strong form of the problem) is a weak solution, but the converse is not always true. While the weak form might look odd, in many mathematical models it could be considered more fundamental than (10.6.2).

10.6.2Galerkin conditions¶

Our goal is to solve a finite-dimensional problem that approximates the weak form of the BVP. Let $\phi_0,\phi_1,\ldots,\phi_m$ be linearly independent functions satisfying $\phi_i(a)=\phi_i(b)=0$ . If we require

\psi(x) = \sum_{i=1}^m z_i \phi_i(x),

(10.6.5)

then (10.6.4) becomes, after some rearrangement,

\begin{split} \sum_{i=1}^m z_i \left[ \int_a^b \bigl[ c(x)u'(x)\phi_i'(x)\,dx + s(x)u(x)\phi_i(x) - f(x) \phi_i(x)\bigr] \, d x \right] = 0. \end{split}

(10.6.6)

One way to satisfy this condition is to ensure that the term inside the brackets is zero for each possible value of $i$ , that is,

\int_a^b \bigl[ c(x)u'(x)\phi_i'(x) + s(x)u(x)\phi_i(x)\bigr] \,dx = \int_a^b f(x)\phi_i(x) \,dx

(10.6.7)

for $i=1,\ldots,m$ . The independence of the $\phi_i$ furthermore guarantees that this is the only possibility, so we no longer need to consider the $z_i$ .

Now that we have approximated the weak form of the BVP by a finite set of constraints, the next step is to represent the approximate solution by a finite set as well. A natural choice is to approximate $u(x)$ the same way as we did the test function $\psi$ , where the $\phi_j$ form a basis for representing the solution:

u(x) = \sum_{j=1}^m w_j \phi_j(x).

(10.6.8)

Substituting (10.6.8) into (10.6.7) implies

\int_a^b \left\{ c(x) \Bigl[ \sum_{j=1}^m w_j \phi_j'(x) \Bigr] \phi_i'(x) + s(x)\Bigl[ \sum_{j=1}^m w_j \phi_j(x) \Bigr]\phi_i(x) \right\} \,dx = \int_a^b f(x)\phi_i(x) \,dx

(10.6.9)

for $i=1,\ldots,m$ . This rearranges easily into

\sum_{j=1}^m w_j \left[ \int_a^b c(x)\phi_i'(x)\phi_j'(x) \,dx + \int_a^b s(x)\phi_i(x)\phi_j(x) \,dx \right] = \int_a^b f(x)\phi_i(x) \,dx,

(10.6.10)

still for each $i=1,\ldots,m$ . These are the Galerkin conditions defining a numerical solution. They follow entirely from the BVP and the choice of the $\phi_i$ .

The conditions (10.6.10) are a linear system of equations for the unknown coefficients $w_j$ . Define $m\times m$ matrices $\mathbf{K}$ and $\mathbf{M}$ , and the vector $\mathbf{f}$ , by

\begin{split} K_{ij} &= \int_a^b c(x)\phi_i'(x)\phi_j'(x) \,dx, \quad i,j=0,\ldots,m,\\ M_{ij} &= \int_a^b s(x)\phi_i(x)\phi_j(x) \,dx, \quad i,j=0,\ldots,m, \\ f_i &= \int_a^b f(x)\phi_i(x) \,dx \quad i=0,\ldots,m. \end{split}

(10.6.11)

Then (10.6.10) is simply

(\mathbf{K}+\mathbf{M})\mathbf{w} = \mathbf{f}.

(10.6.12)

The matrix $\mathbf{K}$ is called the stiffness matrix and $\mathbf{M}$ is called the mass matrix. By their definitions, they are symmetric. The last piece of the puzzle is to make some selection of $\phi_1,\ldots,\phi_m$ and obtain a fully specified algorithm.

Example 10.6.1

Suppose we are given $-u''+4u=x$ , with $u(0)=u(\pi)=0$ . We could choose the basis functions $\phi_k=\sin(kx)$ for $k=1,2,3$ . Then

\begin{split} M_{ij} & = 4 \int_0^\pi \sin(ix) \sin(jx)\, dx, \\ K_{ij} & = ij \int_0^\pi \cos(ix) \cos(jx)\, dx, \\ f_i &= \int_0^\pi x \sin(ix)\, dx. \end{split}

With some calculation, we find

\mathbf{M} = 2\pi \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \qquad \mathbf{K} = \frac{\pi}{2} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 9 \end{bmatrix}, \qquad \mathbf{f} = \pi \begin{bmatrix} 1 \\ -1/2 \\ 1/3 \end{bmatrix}.

Upon solving the resulting diagonal linear system, the approximate solution is

\frac{2}{5}\sin(x) - \frac{1}{8} \sin(2x) + \frac{2}{39}\sin(3x).

10.6.3Finite elements¶

One useful and general choice for the $\phi_i$ are the piecewise linear hat functions constructed in Piecewise linear interpolation. As usual, we select nodes $a=t_0 < t_1 < \cdots < t_n=b$ . Also define

h_i = t_i - t_{i-1}, \qquad i=1,\ldots,n.

(10.6.13)

Then we set $m=n-1$ , and the $\phi_i$ in (10.6.8) are

\phi_i(x) = H_i(x) = \begin{cases} \dfrac{x-t_{i-1}}{h_i} & \text{if $x\in[t_{i-1},t_i]$},\\[2.5ex] \dfrac{t_{i+1}-x}{h_{i+1}} & \text{if $x\in[t_{i},t_{i+1}]$},\\[2.5ex] 0 & \text{otherwise}. \end{cases}

(10.6.14)

Recall that these functions are cardinal, i.e., $H_i(t_i)=1$ and $H_i(t_j)=0$ if $i\neq j$ . Hence

u(x) = \sum_{j=1}^m w_j \phi_j(x) = \sum_{j=1}^{n-1} u_j H_j(x),

(10.6.15)

where, as usual, $u_j$ is the value of the numerical solution at $t_j$ . Note that we omit $H_0$ and $H_n$ , which equal one at the endpoints of the interval, because the boundary conditions on $u$ render them irrelevant.

The importance of the hat function basis in the Galerkin method is that each one is nonzero in only two adjacent intervals. As a result, we shift the focus from integrations over the entire interval in (10.6.11) to integrations over each subinterval, $I_k=[t_{k-1},t_k]$ . Specifically, we use

\begin{split} K_{ij} &= \sum_{k=1}^{n} \left[ \int_{I_k} c(x) H_i'(x) H_j'(x) \,dx\right], \qquad i,j=1,\ldots,n-1, \\ M_{ij} &= \sum_{k=1}^{n} \left[ \int_{I_k} s(x)H_i(x)H_j(x) \,dx\right], \qquad i,j=1,\ldots,n-1, \\ f_i &= \sum_{k=1}^{n} \left[ \int_{I_k} f(x) H_i(x) \,dx\right] \qquad i=1,\ldots,n-1. \end{split}

(10.6.16)

Start with the first subinterval, $I_1$ . The only hat function that is nonzero over $I_1$ is $H_1(x)$ . Thus, the only integrals we need to consider over $I_1$ have $i=j=1$ :

\int_{I_1} c(x) H_1'(x) H_1'(x) \,dx, \qquad \int_{I_1} s(x) H_1(x) H_1(x) \,dx, \qquad \int_{I_1} f(x) H_1(x) \,dx,

(10.6.17)

which contribute to the sums for $K_{11}$ , $M_{11}$ , and $f_1$ , respectively.

Before writing more formulas, we make one more very useful simplification. Unless the coefficient functions $c(x)$ , $s(x)$ , and $f(x)$ specified in the problem are especially simple functions, the natural choice for evaluating all the required integrals is numerical integration, say by the trapezoid formula. As it turns out, though, such integration is not really necessary. The fact that we have approximated the solution of the BVP by a piecewise linear interpolant makes the numerical method second-order accurate overall. It can be proven that the error is still second order if we replace each of the coefficient functions by a constant over $I_k$ , namely the average of the endpoint values:

c(x) \approx \overline{c}_k = \frac{c(t_{k-1})+c(t_k)}{2} \quad \text{for $x\in I_k$}.

(10.6.18)

Thus, the integrals in (10.6.16) can be evaluated solely from the node locations. For instance,

\int_{I_1} c(x) H_1'(x) H_1'(x) \,dx \approx \overline{c}_1 \int_{t_0}^{t_1} h_1^{-2} \, dx = \frac{\overline{c}_1}{h_1}.

(10.6.19)

Now consider interval $I_2=[t_1,t_2]$ . Here both $H_1$ and $H_2$ are nonzero, so there are contributions to all the matrix elements $K_{11}$ , $K_{12}=K_{21}$ , $K_{22}$ , to $M_{11}$ , $M_{12}=M_{21}$ , $M_{22}$ , and to $f_1$ and $f_2$ . Over $I_2$ we have $H_2'= h_2^{-1}$ and $H_{1}'= -h_2^{-1}$ . Hence, the contributions to $K_{11}$ and $K_{22}$ in (10.6.16) are $\overline{c}_2/h_2$ , and the contributions to $K_{12}=K_{21}$ are $-\overline{c}_2/h_2$ . We summarize the relationship by

\frac{\overline{c}_k}{h_k} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \rightsquigarrow \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix},

(10.6.20)

where the squiggly arrow is meant to show that the values of the $2\times 2$ matrix on the left are added to the appropriate submatrix of $\mathbf{K}$ on the right. Similar expressions are obtained for contributions to $\mathbf{M}$ and $\mathbf{f}$ in (10.6.16); see below.

In general, over $I_k$ for $1<k<n$ , we have $H_k'= h_k^{-1}$ and $H_{k-1}'=-h_k^{-1}$ . The stiffness matrix contributions over $I_k$ become

\frac{\overline{c}_k}{h_k} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \rightsquigarrow \begin{bmatrix} K_{k-1,k-1} & K_{k-1,k} \\ K_{k,k-1} & K_{k,k} \end{bmatrix}, \qquad k=2,\ldots,n-1.

(10.6.21)

One finds the contributions to the other structures by similar computations:

\frac{\overline{s}_k h_k}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \rightsquigarrow \begin{bmatrix} M_{k-1,k-1} & M_{k-1,k} \\ M_{k,k-1} & M_{k,k} \end{bmatrix}, \qquad k=2,\ldots,n-1,

(10.6.22)

and

\frac{\overline{f}_k h_k}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \rightsquigarrow \begin{bmatrix} f_{k-1} \\ f_{k} \end{bmatrix}, \qquad k=2,\ldots,n-1.

(10.6.23)

The contribution from $I_n$ affects just $K_{n-1,n-1}$ , $M_{n-1,n-1}$ , and $f_{n-1}$ , and it produces formulas similar to those for $I_1$ .

Each $I_k$ contributes to four elements of each matrix and two of the vector $\mathbf{f}$ , except for $I_1$ and $I_n$ , which each contribute to just one element of each matrix and $\mathbf{f}$ . The spatially localized contributions to the matrices characterize a finite element method (FEM). Putting together all the contributions to (10.6.12) to form the complete algebraic system is often referred to as the assembly process.

10.6.4Implementation and convergence¶

Function 10.6.1 implements the piecewise linear FEM on the linear problem as posed in (10.6.4), using an equispaced grid. The code closely follows the description above.

Algorithm 10.6.1 (fem)

Julia

MATLAB

Python

Piecewise linear finite elements for a linear BVP

fem.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
"""
    fem(c, s, f, a, b, n)

Use a piecewise linear finite element method to solve a two-point
boundary value problem. The ODE is (`c`(x)u')' + `s`(x)u = `f`(x) on
the interval [`a`,`b`], and the boundary values are zero. The
discretization uses `n` equal subintervals.

Return vectors for the nodes and the values of u.
"""
function fem(c, s, f, a, b, n)
    # Define the grid.
    h = (b - a) / n
    x = @. a + h * (0:n)

    # Templates for the subinterval matrix and vector contributions.
    Ke = [1 -1; -1 1]
    Me = (1 / 6) * [2 1; 1 2]
    fe = (1 / 2) * [1; 1]

    # Evaluate coefficent functions and find average values.
    cval = c.(x)
    cbar = (cval[1:n] + cval[2:n+1]) / 2
    sval = s.(x)
    sbar = (sval[1:n] + sval[2:n+1]) / 2
    fval = f.(x)
    fbar = (fval[1:n] + fval[2:n+1]) / 2

    # Assemble global system, one interval at a time.
    K = zeros(n - 1, n - 1)
    M = zeros(n - 1, n - 1)
    f = zeros(n - 1)
    K[1, 1] = cbar[1] / h
    M[1, 1] = sbar[1] * h / 3
    f[1] = fbar[1] * h / 2
    K[n-1, n-1] = cbar[n] / h
    M[n-1, n-1] = sbar[n] * h / 3
    f[n-1] = fbar[n] * h / 2
    for k in 2:n-1
        K[k-1:k, k-1:k] += (cbar[k] / h) * Ke
        M[k-1:k, k-1:k] += (sbar[k] * h) * Me
        f[k-1:k] += (fbar[k] * h) * fe
    end

    # Solve system for the interior values.
    u = (K + M) \ f
    u = [0; u; 0]      # put the boundary values into the result
    return x, u
end

Piecewise linear finite elements for a linear BVP

fem.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function [x, u] = fem(c, s, f, a, b, n)
    %FEM     Piecewise linear finite elements for a linear BVP.
    % Input:
    %   c,s,f    coefficient functions of x describing the ODE (functions) 
    %   a,b      domain of the independent variable (scalars)
    %   n        number of grid subintervals (scalar) 
    % Output:
    %   x        grid points (vector, length n+1)
    %   u        solution values at x (vector, length n+1)

    % Define the grid.
    h = (b - a)/n;
    x = a + h * (0:n)';  

    % Templates for the subinterval matrix and vector contributions.
    Ke = [1 -1; -1 1];
    Me = (1/6) * [2 1; 1 2];
    fe = (1/2) * [1; 1];

    % Evaluate coefficent functions and find average values.
    cval = c(x);   cbar = (cval(1:n) + cval(2:n+1)) / 2;
    sval = s(x);   sbar = (sval(1:n) + sval(2:n+1)) / 2;
    fval = f(x);   fbar = (fval(1:n) + fval(2:n+1)) / 2;

    % Assemble global system, one interval at a time.
    K = zeros(n-1, n-1);  M = zeros(n-1, n-1);  f = zeros(n-1, 1);
    K(1, 1) = cbar(1) / h;
    M(1,1)  = sbar(1) * h / 3;  
    f(1)    = fbar(1) * h / 2;
    K(n-1, n-1) = cbar(n) / h;
    M(n-1, n-1) = sbar(n) * h / 3;
    f(n-1)      = fbar(n) * h / 2;
    for k = 2:n-1
        K(k-1:k, k-1:k) = K(k-1:k, k-1:k) + (cbar(k) / h) * Ke;
        M(k-1:k, k-1:k) = M(k-1:k, k-1:k) + (sbar(k) * h) * Me;
        f(k-1:k)        = f(k-1:k)        + (fbar(k) * h) * fe;
    end  

    % Solve system for the interior values.
    u = (K + M) \ f;
    u = [0; u; 0];      % put the boundary values into the result
end

Piecewise linear finite elements for a linear BVP

fem.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def fem(c, s, f, a, b, n):
    """
    fem(c, s, f, a, b, n)

    Use a piecewise linear finite element method to solve a two-point boundary
    value problem. The ODE is (c(x)u')' + s(x)u = f(x) on the interval
    [a,b], and the boundary values are zero. The discretization uses n equal
    subintervals.

    Return vectors for the nodes and the values of u.
    """
    # Define the grid.
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)

    # Templates for the subinterval matrix and vector contributions.
    Ke = np.array([[1, -1], [-1, 1]])
    Me = (1 / 6) * np.array([[2, 1], [1, 2]])
    fe = (1 / 2) * np.array([1, 1])

    # Evaluate coefficent functions and find average values.
    cval = c(x)
    cbar = (cval[:-1] + cval[1:]) / 2
    sval = s(x)
    sbar = (sval[:-1] + sval[1:]) / 2
    fval = f(x)
    fbar = (fval[:-1] + fval[1:]) / 2

    # Assemble global system, one interval at a time.
    K = np.zeros([n - 1, n - 1])
    M = np.zeros([n - 1, n - 1])
    f = np.zeros(n - 1)
    K[0, 0] = cbar[0] / h
    M[0, 0] = sbar[0] * h / 3
    f[0] = fbar[0] * h / 2
    K[-1, -1] = cbar[-1] / h
    M[-1, -1] = sbar[-1] * h / 3
    f[-1] = fbar[-1] * h / 2
    for k in range(1, n - 1):
        K[k - 1 : k + 1, k - 1 : k + 1] += (cbar[k] / h) * Ke
        M[k - 1 : k + 1, k - 1 : k + 1] += (sbar[k] * h) * Me
        f[k - 1 : k + 1] += (fbar[k] * h) * fe

    # Solve system for the interior values.
    u = np.linalg.solve(K + M, f)
    u = np.hstack([0, u, 0])  # put the boundary values into the result

    return x, u

Example 10.6.2 (Finite element solution of a BVP)

We solve the equation

-(x^2u')' + 4 y = \sin(\pi x), \qquad u(0)=u(1)=0,

in which

c(x) = x^2, \qquad s(x) = 4, \qquad f(x)=\sin(\pi x).

Julia

MATLAB

Python

Example 10.6.2

Here are the coefficient function definitions. Even though $s$ is a constant, it has to be defined as a function for Function 10.6.1 to use it.

c = x -> x^2;
q = x -> 4;
f = x -> sin(π * x);

using Plots
x, u = FNC.fem(c, q, f, 0, 1, 50)
plot(x, u;
    xaxis=(L"x"),  yaxis = (L"u"),
    title = "Solution by finite elements", legend=:none)

Example 10.6.2

Here are the coefficient function definitions. Even though $s$ is a constant, it has to be defined as a function for Function 10.6.1 to use it.

c = @(x) x.^2;
q = @(x) 4 * ones(size(x));
f = @(x) sin(pi * x);

[x,u] = fem(c, q, f, 0, 1, 50);
clf,  plot(x, u)
xlabel('x'),  ylabel('u')
title('Solution by finite elements')

Example 10.6.2

Here are the coefficient function definitions. Even though $s$ is a constant, it has to be defined as a function for Function 10.6.1 to use it.

c = lambda x: x**2
q = lambda x: 4 * ones(len(x))
f = lambda x: sin(pi * x)

x, u = FNC.fem(c, q, f, 0, 1, 50)
plot(x, u)
xlabel("$x$"),  ylabel("$u$")
title("Solution by finite elements");

Because piecewise linear interpolation on a uniform grid of size $h$ is $O(h^2)$ accurate, the accuracy of the FEM method based on linear interpolation as implemented here is similar to the second-order finite-difference method.

10.6.5Exercises¶

Exercise 10.6.1

⌨ For each linear BVP, use Function 10.6.1 to solve the problem and plot the solution for $n=40$ . Then for each $n=10,20,40,\ldots,640$ , compute the norm of the error. Make a log-log convergence plot of error versus $n$ and compare graphically to second-order convergence.

(a) $-u''+u=-8 + 16 x^2 - x^4, \quad u(0) =u(2) =0$

Exact solution: $x^2(4-x^2)$

(b) $[(2+x)u']' +11x u = -e^x \left(12 x^3+7 x^2+1\right), \quad u(-1) =u(1) =0$

Exact solution: $e^x \left(1-x^2\right)$

(c) $u''+x(u'+u) = -x[4 \sin(x)+5 x \cos(x)], \quad u(0) =u(2\pi) =0$

Exact solution: $-x^2\sin(x)$

Exercise 10.6.5

Suppose the Dirichlet boundary conditions $u(a)=u(b)=0$ are replaced by the homogeneous Neumann conditions $u'(a)=u'(b)=0$ .

(a) ✍ Explain why the weak form (10.6.4) can be derived without any boundary conditions on the test function $\psi$ .

(b) ⌨ The result of part (a) suggests replacing (10.6.15) with

u(x) = \sum_{j=0}^{n} u_j H_j(x)

(10.6.24)

and making (10.6.16) hold for all $i,j$ from 0 to $n$ . Modify Function 10.6.1 to do this and thereby solve the Neumann problem. (Note that $I_1$ and $I_n$ now each make multiple contributions, like all the other integration subintervals.)

(c) ⌨ Test your function on the problem

u''+u = -2 \sin(x), \quad u'(0)=u'(1)=0,

(10.6.25)

whose exact solution is $(x-1)\cos(x) - \sin(x)$ . Show second-order convergence.

Preface

Nonlinearity and boundary conditions

Preface

Next steps