06. Discretisation

In one spatial dimension $x$, consider the boundary value problem (BVP) with solution $u(x)$:

$$A(x) \frac{\partial^2 u}{\partial x^2}+B(x) \frac{\partial u}{\partial x}+C(x) u = D(x),$$

where $A(x), B(x), C(x)$ and $D(x)$ are functions. This equation is supplemented with Dirichlet or Neumann boundary conditions.

We will solve this problem by discretising the equation using finite differences to turn the differential equation into an algebraic one with matrices.

This will involve inverting a matrix which we will achieve with the linear algebra package LAPACK.

Solving such equations are described in this exemplar, but we will list the important steps here:

1. Preprocessing:
   • Using reader.f90 we read in the specified input data from the file settings.input.
   • Using domain.f90 we build the uniform computational grids for the problem.
   • With equations/definitions.f90 we then specify the equations we are solving, with boundary conditions.
2. Discretisation:
   • In equations.f90, specifically in the subroutine build_the_matrix, the equation is spatially discretised into discrete algebraic equations using finite differences.
   • Both the left and right hand side of the equation above need to be assembled, including all boundary conditions.
   • The equation is then transformed in a banded form. This is due to the fact that most of the matrix entries will be zero when finite differences are used. Again refer to this exemplar for more details.
3. Solution: we then solve the problem in the module solve_bvp.f90 by calling on the LAPACK linear solver defined in linear_algebra.f90 with the subroutine solver_banded_double_precision.
4. The solution is then outputted to a data file called BVP.dat.