Computational Linear Algebra
Python code for finite elements polynomial approximation, ODE and PDE solvers
––– views

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- 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 , where is time and is space. This equation can be solved quickly by transforming this second order PDE to a system of first order PDEs
and then discretising using the implicit midpoint rule. The resulting matrix-vector system can be written as , where 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
where is tridiagonal, for some vectors to be determined. Using repeated applications of the Sherman-Morrison formula, can be written in terms of , and the resulting expression for (the solution to the equation) can be simplified considerably using the known forms of . An efficient algorithm can therefore construct an solution to this PDE, compared to from LU factorisation.
The code in this folder solves the system 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 for each coordinate. This is done by:
- Constructing a new matrix-vector system to solve for all timesteps and positions at once
- Writing for banded matrices , where 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