Week 19: Finite difference methods

Overview

This week will introduce the idea of finite difference methods and how they can be used to solve boundary value problems (BVPs). We will first revisit how to perform numerical differentiation on a computer using finite differences. We will then see how the finite difference method can be used to solve BVPs for linear and nonlinear ordinary differential equations (ODEs). The code that you develop this week will form an essential role in the upcoming weeks when the focus shifts towards numerically solving partial differential equations.

Supplementary material

Use the links below to find additional notes on

• Finite-difference formulae: derivation and error analysis

• Solving the nonlinear Poisson equation with Newton's method

Exercise

The goal of the exercise this week is to create a BVP solver that is capable of finding numerical solutions to ODEs of the form

where the domain of the problem is given by $a \leq x \leq b$. In this equation, $D > 0$ is a parameter. The term $q(x, u; \mu)$ represents a function that depends on the spatial coordinate $x$, the solution $u(x)$, and a parameter $\mu$. Your solver should be able to handle all three types of boundary conditions (Dirichlet, Neumann, and Robin).

You should design your BVP solver so that it is compatible with the numerical continuation code you developed in Week 17. This will allow you to use numerical continuation to track how the solution $u(x)$ changes under the variation of one of the parameters in the problem.

Ideally, the user should be able to select how the algebraic system that arises from discretising the problem will be solved. For example, the discrete problem could be solved using SciPy's root function, NumPy's solve function (if the problem is linear), or your own implementation of Newton's method (see the supplementary notes above on how to implement this).

Steps

There is a lot going on in this problem. These steps are designed so that you can see how to build your code up from simpler problems. Test your code extensively and only add complexity once you are certain it works correctly.

1. Use finite differences to find a numerical solution to

   $$\frac{\mathrm{d}^2 u}{\mathrm{d} x^2} = 0, \quad u(a) = \gamma_1, \quad u(b) = \gamma_2.$$

   In this case, the exact solution to the problem is given by

   $$u(x) = \left(\frac{\gamma_2 - \gamma_1}{b-a}\right)(x - a) + \gamma_1,$$

   which you can use to test your code.

2. Extend your code so that it can account for a source term in the ODE:

   $$D \frac{\mathrm{d}^2 u}{\mathrm{d} x^2} + q(x) = 0, \quad u(a) = \gamma_1, \quad u(b) = \gamma_2.$$

   Hint: The simplest place to start is to set $q(x) = 1$; in this case, the exact solution is given by

   $$u(x) = -\frac{1}{2D}(x-a)(x-b) + \left(\frac{\gamma_2 - \gamma_1}{b-a}\right)(x - a) + \gamma_1.$$

   Once you've developed code for the problem with $q(x) = 1$, add an $x$ dependence into $q$. Can you find some exact solutions for this case?

3. Generalise your code so that the source term $q$ can now depend on the solution $u$ as well as a parameter $\mu$. Use your code to solve the Bratu problem

   $$\begin{aligned} D \frac{\mathrm{d}^2 u}{\mathrm{d} x^2} + e^{\mu u} = 0, \qquad u(0) = u(1) = 0, \end{aligned}$$

   when $D = 1.0$ and $\mu = 0.1$. Plot the solution $u$ as a function of $x$. The Bratu problem appears in mathematical models of combustion and thermal runaway, in which case $u$ is the temperature.

   Hint: If $q$ depends nonlinearly on the solution $u$, then a good initial guess of the solution is usually required for the nonlinear solver (e.g. SciPy's root function or Newton's method) to converge. For this problem, a good initial guess can be found by noting that when $\mu$ is small, the exponential can be approximated as $e^{\mu u} \approx 1$. The solution in this case is given by (5), which can be used to form an initial guess.

4. Update your code so that it can account for Dirichlet, Neumann, or Robin boundary conditions.