Project 2: Iterative solvers, FEM, heat equation#

Topics covered: Iterative solvers, finite element methods.

Programming instructions#

You are welcome to use any of the codes that have been provided as part of the course.

Report instructions#

Your written report needs to include the derivations and plots that are requested under “Assignment 1-6” below. There are also a few questions that you need to answer.
Submit your report as one pdf that includes the source code.
In addition to the pdf, submit your source code files separately so that we can download them and test your code.
All equations must be nicely formatted and all figures must have captions. The text must refer to the figures and further explain what they illustrate.
The report does not need to have a specific structure (Introduction, Methods, Results, etc.). You are free to select your own headings as long as you include the requested derivations, plots, answers, and code.
We recommend using LaTeX to generate the report. There are tutorials on e.g. Overleaf. It is also possible (but more cumbersome) to generate a nicely formatted report using e.g. Word.
Put all group members’ names on the first page.

Solving linear systems#

Consider a sequence of \(N\) linear systems,

\[ Ax_i = b_i, \quad i = 1, 2, \ldots, N \]

where \(A\) is a known, sparse \(n \times n\) matrix, the \(b_i\) are known \(n \times 1\) vectors, and the \(x_i\) are unknown \(n \times 1\) solution vectors. You will implement and evaluate four different ways to solve this sequence of systems:

Gaussian elimination, using for example scipy.sparse.linalg.spsolve() or MATLAB’s backslash operator.
LU factorization
The Jacobi method
The conjugate gradient method

The first step is to implement methods 1-4 and verify that they work as expected (we will assume that method 1 prdouces the correct answer). For this test, we consider the system matrix

\[\begin{split} A = \begin{bmatrix} 4 & 1 & & & & \ldots \\ 1 & 4 & 1 & & & \ldots \\ & 1 & 4 & 1 & & \ldots \\ & & \ddots & \ddots & \ddots & \\ & & & 1 & 4 & 1 \\ & & & & 1 & 4 \\ \end{bmatrix} \end{split}\]

A Python function that produces this matrix is provided below.

import numpy as np
import scipy.sparse as spsp
def system_matrix(n):
    A = 1*np.diag(np.ones(n-1), 1) + 1 * \
        np.diag(np.ones(n-1), -1) + 4*np.diag(np.ones(n), 0)
    A = spsp.csr_matrix(A)
    return A

LU factorization and the Jacobi method perform some of the computations “once and for all” in a pre-processing step. Those computations are then used to compute the \(x_i\) more efficiently. For the Jacobi method, you should finish the implementation of the following two functions:

def jacobi_setup(A):
   """Splits the matrix A and returns appropriate matrices"""

and

def jacobi_solve(Dinv, L_plus_U, b, x0=None, tol=1e-6):
   """Solves a linear system Ax=b using the Jacobi method"""

   # Dinv:     inv(D), where D is the diagonal part of A
   # L_plus_U: L+U (upper and lower triangular parts of A)
   # b:        right-hand side vector
   # x0:       Initial guess (if None, the zero vector is used)
   # tol:      Relative error tolerance

   # If no initial guess supplied, use the zero vector
   if x0 is None:
       x_new = 0*b
   else:
       x_new = x0

   # More code here ...

   # Return approximate solution and number of iterations required
   return x, n_iter

To check whether approximation \(x_k\) is accurate enough, you may approximate the relative error according to

\[ e_k \approx \frac{\lVert x_k - x_{k-1} \rVert_2}{\lVert x_k \rVert_2} \]

Assignment 1#

Complete the functions jacobi_setup and jacobi_solve. Then use them to run the function test_solver, provided below. The function solves \(N\) \(n \times n\) systems with random right-hand sides and compares with the solution obtained by Gaussian elimination. If the error is larger than \(10*tol\), an error message is printed.

You may use the following implementation (or a ported version of it) of the conjugate gradient method:

After completing jacobi_setup and jacobi_solve, run test_solver with N=3, n=1000, tol=1e-6. Run with method set to "lu", "cg", and "jacobi", to verify that all solvers are working as expected.

NOTE: splu and spsolve may throw warnings related to the format of A and/or b. You may ignore these warnings.

You do NOT need to report anything from this assignment, but it is important to verify that the solvers are working as expected as you will need them in the next steps.

Assignment 2#

Run the function run_and_time, provided below, with n=1e4, N=100, tol=1e-6 and methods = ['cg', 'gauss', 'lu', 'jacobi'].

The function will print the time required for the solution process, disregarding setup time. Report the results in a table of this form:

Solver	Time (s)	Total iterations
Gaussian elimination	x	N/A
LU factorization	x	N/A
Jacobi	x	x
Conjugate gradient	x	x

Which is the fastest direct method? Which is the fastest iterative method? Which is the fastest out of all of them? Briefly discuss the results.

A FEM solver for the heat equation#

Consider the heat equation with constant coefficients and homogeneous Dirichlet boundary conditions:

\[\begin{split} \begin{array}{rll} u_t = au_{xx}, & x \in (0, 1), & 0<t<T,\\ u=0, & x=0, & 0<t<T, \\ u=0, & x=1, & 0<t<T, \\ u = f, & x \in (0, 1), & t=0,\\ \end{array} \quad ({\color{red}*}) \end{split}\]

where \(a>0\) is a real constant.

Assignment 3#

Derive the weak form of \(({\color{red}*})\) with appropriate spaces.

Assignment 4#

Derive a finite element approximation of the weak form, using appropriate spaces of piecewise linear functions on a uniform mesh of \(n\) intervals of length \(h=1/n\). That is, the nodes are

\[ x_j = jh, \quad j=0,1,\ldots,n. \]

Assignment 5#

Derive the system of ODEs that the finite elment discretization results in. Evaluate all integrals in the mass and stiffness matrices exactly.

Assignment 6#

Implement the finite element method derived in Assignment 5, combined with the 4th order Runge–Kutta method for time integration. Note that you will have to solve a linear system involving the mass matrix in every Runge–Kutta stage. Implement two ways of solving the system: the fastest direct method and the fastest iterative method, based on your results in Assignment 2.

To impose the initial condition \(u(x,0) = f(x)\), you may impose that the finite element solution \(u_h\) should equal \(f\) at the nodes (there are other ways to impose the initial condition, but this is arguably the simplest). If the discrete solution is

\[ u_h(x,t) = \sum \limits_{j=1}^{n-1} \xi_j(t) \varphi_j(x), \]

where the \(\varphi_j\) are the usual hat functions, then the initial condition becomes

\[ \xi_j(0) = f(x_j), \quad j = 1,2,\ldots,n-1. \]

Consider the following rapidly oscillating initial data:

\[ f(x) = \sin(kx), \]

with \(k=100 \cdot 2\pi\). The corresponding exact solution is

\[ u(x,t) = \sin(kx)e^{-ak^2t}. \]

Set \(a=1\) and the final time \(T=10^{-5}\). Select the time step according to

\[ \Delta t = 0.1 \sqrt{a} h^2. \]

For your iterative solver, use the relative error tolerance tol=1e-6. The error at node \(j\) is

\[ e_j = \xi_j(T) - u(x_j, T). \]

Use the 2-norm to measure the error, so that the relative error is

\[ e = \sqrt{\frac{\sum \limits_{j=1}^{n-1} e_j^2}{\sum \limits_{j=1}^{n-1} u(x_j, T)^2}}. \]

Run simulations with your two solvers, using \(n=2000\), \(n=4000\), and \(n=8000\). Present the results in two tables like the one below, one for each solver.

n	Time (s)	Relative error
2000	x	x
4000	x	x
8000	x	x

Measure the execution time for the time-stepping loop only. Disregard setup time and switch of plotting and other non-essential operations.

Does the error converge with the expected rate?

Scientific Computing for PDE

Project 2: Iterative solvers, FEM, heat equation

Contents