Computational Linear Algebra

Python code for finite elements polynomial approximation, ODE and PDE solvers

––– views

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-pp 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 utt=uxxu_{tt} = u_{xx}, where tt is time and xx is space. This equation can be solved quickly by transforming this second order PDE to a system of first order PDEs

wtuxx=0 and utw=0w_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=bAx = b, where AA 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+u1v1T+u2v2TA = T + u_1 v_1^T + u_2 v_2^T

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

The code in this folder solves the system Ax=bAx=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 ui,wiu_i, w_i for each coordinate. This is done by:

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