Computational Linear Algebra

This post contains a selection of work completed for the 2020/2021 Computational Linear Algebra course at Imperial. The course focused on numerical methods to solve linear systems in Python, including:

QR and LU decomposition
analysis and modifications of algorithm to utilise matrix structures
eigenvalue problems
iterative methods

My work on this course can be found here. A selection of particularly interesting problems are discussed below.

Piecewise polynomial approximation for the pressure-depth relationship of porous rock

Consider a porous rock below the ground, such that the geology of the rock changes below depth 1 (in non-dimensional units). The relationship between depth and pressure will therefore also change at depth 1, but physical laws tell us that pressure should be continuous across the boundary. Making an assumption on the change in derivative at the boundary allows construction of a degree- $p$ polynomial separately for the rock above and below depth 1. This can be formulated as a least squares problem, and, with a suitable change of variables, can be solved with QR factorisation. A similar approach can be used to approximate smooth, closed curves from the unit circle to the 2D plane.

Solutions to the 1-dimensional wave equation with periodic boundary conditions

Individual timesteps

It can be shown that 1-dimensional waves (e.g. a pulse travelling along a string) satisfy the partial differential equation $u_{tt} = u_{xx}$ , where $t$ is time and $x$ is space. This equation can be solved quickly by transforming this second order PDE to a system of first order PDEs

w_t - u_{xx} = 0 \text{ and } u_t - w = 0

and then discretising using the implicit midpoint rule. The resulting matrix-vector system can be written as $Ax = b$ , where $A$ is tridiagonal with non-zero entries in the top right and bottom left corners. As a result, banded matrix algorithms are not helpful in this case, but instead we can write

A = T + u_1 v_1^T + u_2 v_2^T

where $T$ is tridiagonal, for some vectors $u_i, v_i, i = 1, 2$ to be determined. Using repeated applications of the Sherman-Morrison formula, $A^{-1}$ can be written in terms of $T^{-1}$ , and the resulting expression for $x$ (the solution to the equation) can be simplified considerably using the known forms of $u_i, v_i$ . An efficient algorithm can therefore construct an $O(m)$ solution to this PDE, compared to $O(m^3)$ from LU factorisation.

The code in this folder solves the system $Ax=b$ efficiently, plots the solution and provides tests to demonstrate the algorithm works correctly.

Parallel timesteps

The approach above can be modified to simultaneously solve for each timestep of $u_i, w_i$ for each coordinate. This is done by:

Constructing a new matrix-vector system $AU = r$ to solve for all timesteps and positions at once
Writing $A = C_1 \otimes I + C_2 \otimes B$ for banded matrices $C_1, C_2$ , where $\otimes$ is the tensor product
Identifying a change of variables to make the system block diagonal
Using discrete Fourier (and inverse Fourier) transforms to apply an iterative scheme which converges rapidly to the solution