Lab 2: Finite element methods#

Topics covered: Weak form, spaces of piecewise polynomial functions.

Note: The labs are not mandatory, but highly recommended.

Introduction#

The finite element method is a general technique for computing numerical solutions to differential equations. In this computer session we will study the derivation and implementation of a finite element method for a two-point boundary value problem, namely,

\[\begin{split} \begin{aligned} -u''&=f, \quad a<x<b, \quad ({\color{red}1}) \\ u(a)=u(b)&=0, \end{aligned} \end{split}\]

where \(u\) is the sought solution and \(f\) a given function.

The weak form#

A finite element method is always derived from the weak or variational formulation of the problem at hand. Let us define the vector space

\[ H^1_0=\{ v: \| v \|^2 + \| v' \|^2 < \infty, ~ v(a)=v(b)=0\}. \]

Multiplying \(-u''=f\) by a function \(v\in H^1_0\) and integrating, using the integration-by-parts formula, gives

\[\begin{split} \begin{aligned} \int_a^b fv\,dx &= \int_a^b -u''v \\ &= -[u'v]_a^b + \int_a^b u'v'\,dx \\ &= -u'(b)v(b) + u'(a)v(a) + \int_a^b u'v'\,dx. \end{aligned} \end{split}\]

Since \(v(a)=v(b)=0\) the boundary terms vanish and we are left with

\[ \begin{aligned} \int_a^b fv\,dx = \int_a^b u'v' \,dx. \end{aligned} \]

Hence, the weak or variational form of \(({\color{red}1})\) reads:

Find \(u\in H^1_0\), such that

\[ \begin{aligned} \label{VF} \int_a^b u' v' \,dx = \int_a^b fv\,dx, \quad ({\color{red}2}) \end{aligned} \]

for all \(v\in H^1_0\).

Finite element approximation#

The space \(H^1_0\) contains many functions and it is therefore just as hard to find a function \(u\in H^1_0\) which satisfies the variational equation \(({\color{red}2})\) as it is to solve the original problem \(({\color{red}1})\). The approximation of the finite element method is therefore to look for a solution approximation \(u_h\) within a small (finite dimensional) subspace \(V_h\) of \(H^1_0\), typically a space of piecewise polynomials.

Let \(V_h\) be the space of all continuous piecewise linear functions, which vanish at the end points \(a\) and \(b\), on a uniform partition \(a = x_0 < x_1 < x_2 < \ldots < x_N = b\) of the interval \(a\leq x\leq b\) into \(N\) subintervals of equal length \(h\). Moreover let \(\{ \varphi_j\}\) be the set of hat basis functions of \(V_h\) associated with the \(N-1\) nodes \(x_j\), \(j = 1,2\ldots, N-1\), such that

\[\begin{split} \begin{aligned} \varphi_i(x_j) = \left \{ \begin{array}{l} 1, \quad \mbox{if } i = j, \\ 0, \quad \mbox{if } i \neq j. \end{array} \right . \end{aligned} \end{split}\]

Note that there are no hat functions \(\varphi_0\) and \(\varphi_N\) since the functions of \(V_h\) must vanish at the end points \(x_0\) and \(x_N\). Otherwise \(V_h\) would not be a subspace of \(H^1_0\).

The finite element approximation of \(({\color{red}2})\) thus reads:

Find \(u_h \in V_h\) such that

\[ \begin{aligned} \label{FEM} \int_a^b u_h' v' \,dx = \int_a^b fv\,dx, \quad ({\color{red}3}) \end{aligned} \]

for all \(v\in V_h\).

It can be shown that \(({\color{red}3})\) is equivalent to the \(N-1\) equations

\[ \begin{aligned} \int_a^b u_h' \varphi_i' \,dx = \int_a^b f\varphi_i\,dx, \quad i=1,2,\ldots, N-1, \end{aligned} \]

because if \(({\color{red}3})\) is satisfied for each basis function \(\varphi_j\) separately, then it is also satisfied for any linear combination of them, and hence any function \(v \in V_h\).

Expanding \(u_h\) as a linear combination of hat functions we make the ansatz

\[ \begin{aligned} \label{u_h} u_h = \sum_{j=1}^{N-1} \xi_j \varphi_j(x), \quad ({\color{red}4}) \end{aligned} \]

where \(\xi_j\), \(j=1,2,\ldots,N-1\) are \(N-1\) coefficients, the nodal values of \(u_h\), to be determined.

Inserting \(({\color{red}4})\) in \(({\color{red}3})\) yields

\[ \begin{aligned} \sum_{j=1}^{N-1} \xi_j \int_a^b \varphi_j' \varphi_i'\,dx = \int_a^b f \varphi_i \,dx \quad i=1,2,\ldots,N-1, \end{aligned} \]

which is an \((N-1) \times (N-1)\) system of equations for \(\xi_j\). In matrix form we write

\[ \begin{aligned} \label{Ax=b} A \xi = b, \quad (\color{red}5) \end{aligned} \]

where \(A\) is an \((N-1)\times (N-1)\) matrix, the so-called stiffness matrix, with entries

\[ \label{stiffnessmatrix} A_{i,j} = \int^b_a \varphi'_i(x) \, \varphi'_j(x) \, dx, \quad i,j=1,2,\ldots,N-1, \]

\(\xi=(\xi_1,\xi_2,\ldots,\xi_{N-1})^T\) is an \((N-1)\) vector containing the unknown coefficients \(\xi_j\), \(j=1,2,\ldots,N-1\), and \(b\) is an \((N-1)\) vector, the so-called load vector, with entries

\[ \label{loadvector} b_i=\int_a^b f(x)\varphi_i(x)\,dx, \quad i=1,2,\ldots,N-1. \]

Computer Implementation#

Assembling the stiffness matrix and load vector#

The explicit expression for a hat function \(\varphi_i(x)\) is given by

\[\begin{split} \begin{aligned} \label{hatfcn} \varphi_i(x) = \left\{ \begin{array}{ll} (x-x_{i-1})/h, &\mbox{if } x_{i-1} \leq x \leq x_i,\\ (x_{i+1}-x)/h, & \mbox{if } x_i \leq x \leq x_{i+1}, \quad ({\color{red}6})\\ 0, & \mbox{otherwise.} \end{array} \right. \end{aligned} \end{split}\]

Hence the derivative \(\varphi'_i\) takes the values \(+1/h\), \(-1/h\), or \(0\), depending on the interval.

Using \(({\color{red}6})\) we can calculate the entries of the stiffness matrix. For \(|i - j| > 1\), we have \(A_{i,j}=0\), since \(\varphi_i\) and \(\varphi_j\) lack overlapping support. However, if \(i = j\), then

\[\begin{split} \begin{aligned} \int_a^b \varphi_i'^2\,dx &= \int_{x_{i-1}}^{x_{i}} \varphi_i'^2 \,dx + \int_{x_{i}}^{x_{i+1}} \varphi_i'^2 \,dx \\ &= \int_{x_{i}}^{x_{i+1}} \frac{1}{h^2} \,dx + \int_{x_{i}}^{x_{i+1}} \frac{(-1)^2}{h^2} \,dx = \frac{2}{h}. \end{aligned} \end{split}\]

where we have used that \(x_i-x_{i-1}=x_{i+1}-x_i=h\). Further, if \(j = i + 1\), then

\[\begin{split} \begin{aligned} A_{i,i+1} &= \int_a^b \varphi_i'\varphi_{i+1}'\,dx = \int_{x_{i}}^{x_{i+1}} \varphi_i'\varphi_{i+1}'\,dx \\ &= \int_{x_{i}}^{x_{i+1}} \frac{(-1)}{h} \cdot \frac{1}{h}\,dx = -\frac{1}{h}. \end{aligned} \end{split}\]

Similarly, we also have \(A_{i,i-1}=-1/h\). Thus the stiffness matrix is:

\[\begin{split} \begin{aligned}\label{original} A=\frac{1}{h} \begin{bmatrix} 2 & -1 & & & &\\ -1 & 2 & -1 & & &\\ & -1 & 2 & -1 & &\\ & & \ddots & \ddots & \ddots &\\ & & & -1 & 2 & -1\\ & & & & -1 & 2 \end{bmatrix} .\quad ({\color{red} 7}) \end{aligned} \end{split}\]

The entries \(b_i\) of the load vector must often be evaluated using quadrature, since they involve the function \(f\) which can be hard to integrate analytically. For example, using the trapezoidal rule one obtains the approximate load vector entries

\[ \begin{aligned} \label{trapzload} b_i = \int_a^b f \varphi_i\,dx = \int_{x_{i-1}}^{x_{i+1}} f \varphi_i \,dx \approx f(x_i)h. \end{aligned} \]

Assembly#

It is sometimes convenient (but perhaps not advisable) to rewrite the linear system of equations \((\color{red}5)\) as

\[\begin{split} \begin{aligned} \label{extended} \frac{1}{h} \begin{bmatrix} 10^6 & -1 & & & &\\ -1 & 2 & -1 & & &\\ & -1 & 2 & -1 & &\\ & & \ddots & \ddots & \ddots &\\ & & & -1 & 2 & -1\\ & & & & -1 & 10^6 \end{bmatrix} \begin{bmatrix} \xi_0\\ \xi_1\\ \xi_2\\ \vdots\\ \xi_{N-1}\\ \xi_N \end{bmatrix} = h \begin{bmatrix} f(x_0)\\ f(x_1)\\ f(x_2)\\ \vdots\\ f(x_{N-1})\\ f(x_N) \end{bmatrix} , \end{aligned} \end{split}\]

that is, to extend the original \((N-1)\times (N-1)\) system of equations into an \((N+1) \times (N+1)\) system of equations by adding additional equations for the nodal values \(\xi_0=u_h(x_0)\) and \(\xi_N=u_h(x_N)\). Recalling the boundary conditions \(u(a)=u(b)=0\) it is clear that we want \(\xi_0=\xi_N=0\). This is accomplished by setting the diagonal entries \(A_{0,0}\) and \(A_{N+1,N+1}\) to a big number, say \(10^6\). The entries \(b_0\) and \(b_N\) of the extended load vector are superfluous and can be chosen arbitrarily.

The extended system of equations can be computed using the standard assembly technique. That is, the stiffness matrix is successively computed by looping over the elements and adding each element’s contribution to the global stiffness matrix:

\[\begin{split} \begin{aligned} \nonumber \frac{1}{h} \begin{bmatrix} 1 & -1 &~ & & & \\ -1 & 1 &~ & & & \\ ~&~ &~ & & & \\ ~&~ &~ & & & \\ ~&~ &~ & & & \\ ~&~ &~ & & & \end{bmatrix} + \frac{1}{h} \begin{bmatrix} ~ & ~ & ~ & & & \\ ~ & 1 & -1 & & & \\ ~ & -1 & 1 & & & \\ ~ &~ &~ & & & \\ ~ &~ &~ & & & \\ ~ &~ &~ & & & \end{bmatrix} + \ldots + \frac{1}{h} \begin{bmatrix} & & ~ & & ~ & ~ \\ & & ~ & & ~ & ~ \\ & & ~ & & ~ & ~ \\ & & ~ & & ~ & ~ \\ & & ~ & & 1 & -1 \\ & & ~ & & -1 & 1 \\ \end{bmatrix} . \end{aligned} \end{split}\]

Each element matrix is computed by integrating over a single element.

A routine for assembling the stiffness matrix is listed below.

def my_stiffness_matrix_assembler(x):
    #
    # Returns the assembled stiffness matrix A.
    # Input is a vector x of node coords.
    #
    N = len(x) - 1                  # number of elements
    A = spsp.dok_matrix((N+1, N+1)) # initialize stiffness matrix
    for i in range(N):              # loop over elements
        h = x[i+1] - x[i]           # element length
        A[i, i] += 1/h              # assemble element stiffness
        A[i, i+1] += -1/h
        A[i+1, i] += -1/h
        A[i+1, i+1] += 1/h
    A[0, 0] = 1.e+6                 # adjust for BC
    A[N, N] = 1.e+6
    return A.tocsr()

The load vector is computed analogously.

def my_load_vector_assembler(x):
    #
    # Returns the assembled load vector b.
    # Input is a vector x of node coords.
    #
    N = len(x) - 1
    B = np.zeros(N+1)
    for i in range(N):
        h = x[i+1] - x[i]
        B[i] = B[i] + f(x[i])*h/2
        B[i+1] = B[i+1] + f(x[i+1])*h/2
    return B

def f(x):
    return 2

Note that these algorithms also work when the node points are not uniformly distributed.

Putting it all together we get the main routine:

def main():
    a = 0                                 # left end point of interval
    b = 1                                 # right
    N = 5                                 # number of intervals
    h = (b-a)/N                           # mesh size
    x = np.arange(a, b+h/2, h)            # node coords
    A = my_stiffness_matrix_assembler(x)
    B = my_load_vector_assembler(x)
    xi = splg.spsolve(A, B)               # solve system of equations
    plt.plot(x, xi)                       # plot solution
    plt.xlabel('x')
    plt.show()

Exercise 1: Implement the finite element solver outlined above. A Python implementation is available here: Lab2.py. Test the code by solving

\[\begin{split} \begin{aligned} -u''&=2, \quad 0<x<1, \\ u(0)=u(1)&=0. \end{aligned} \end{split}\]

Use uniform meshes with \(h=1/2\), \(1/4\), \(1/16\) and \(1/256\). Compare with the exact solution \(u=x(1-x)\).

Exercise 2: Modify your code so that it can handle the inhomogeneous boundary condition \(u(1)=g\). Test your code by solving

\[\begin{split} \begin{aligned} -u''&=0, \quad 0<x<1, \\ u(0)&=0, \\ u(1)&=7. \end{aligned} \end{split}\]

Hint: Can you determine the exact solution? It is a simple function!

Exercise 3: Investigate the convergence rate of the finite element approximation \(u_h\). Set up a problem for which the exact solution \(u\) is known and study how the error \(e=u-u_h\) behaves when the mesh size \(h\) decreases. Use the so-called energy norm \(\|\cdot\|_E\) to measure the error \(e\). The energy norm depends on the problem, but for the model problem \(({\color{red} 1})\) it is defined by

\[ \|w\|_E^2=\int_a^b (w')^2\,dx. \]

To compute \(\|e\|_E\), note that \(\|e\|_E^2=\|u'\|^2-\|u_h'\|^2\) due to orthogonality. Further, the term \(\|u_h'\|^2\) may be numerically evaluated as \(\xi^T A \xi\).

Note: Do not use a linear function as the exact solution, since \(u_h\) will be exact.

Problem 4: In the extended stiffness matrix, we added large values (\(10^6\)) in the top-left and bottom-right corners. This is similar to enforcing a Robin BC:

\[ a u' = \kappa(u-g) , \]

and using a large value of \(\kappa\) so that the Robin BC approximates a Dirichlet BC. This approach causes unnecessary stiffness, and fortunately there is a better way to enforce Dirichlet BC.

Assume that the boundary conditions are given as \(u(a) = g_1\) and \(u(b) = g_2\), where \(g_1\) and \(g_2\) are given boundary values. We may replace the first and last rows of the matrix \(A\) by the corresponding rows of the identity matrix, and replace the first and last entries of the load vector by the corresponding boundary values, i.e. \(b_0=g_1\) and \(b_N=g_2\). Solve this system and compare this solution to the one you got in Problem 1. Use the same right-hand-side \(f=2\). Compare the condition numbers (np.linalg.cond(A) in Python, cond(A) in Matlab) of the stiffness matrices in the two computations. We will see later that a large condition number is a bad thing when solving linear systems with iterative methods. It is important to discretize the PDE in a way that leads to a small condition number.

Scientific Computing for PDE

Lab 2: Finite element methods

Contents