Week 2: ODE BVPs

\[ \renewcommand{\vec}[1]{\boldsymbol{#1}} \newcommand{\td}[2]{\frac{\mathrm{d}#1}{\mathrm{d}#2}} \newcommand{\tdd}[2]{\frac{\mathrm{d}^2#1}{\mathrm{d}#2^2}} \newcommand{\pd}[2]{\frac{\partial#1}{\partial#2}} \newcommand{\pdd}[2]{\frac{\partial^2#1}{\partial#2^2}} \]

Overview

This week showcase how finite difference methods can be used to solve boundary value problems (BVPs) for ODEs. In short, the finite difference method converts ODEs into algebraic systems of equations, which can then be solved using two main approaches:

  1. Writing your own solvers using linear algebra functions.
  2. Using built-in solvers provided by Python (SciPy) and Matlab.

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

Exercises

  1. Use finite differences to find a numerical solution to \[ \tdd{u}{x} = 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, \tag{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 \tdd{u}{x} + 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 that can be used for code validation?

  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 equation \[ \tdd{\phi}{x} + \mu(\phi - \phi^3) = 0 \] with the boundary conditions \(\phi(-1) = -1\) and \(\phi(1) = 1\) when (a) \(\mu = 0.5\) and (b) \(\mu = 100\).

    Hint 1: If \(q\) depends nonlinearly on the solution, as with this problem, then a good initial guess of the solution is usually required for the nonlinear solver to converge. For this problem, a good initial guess can be found by noting that when \(\mu\) is small, the source term \(q\) is approximately zero. The solution in this case is given by (1), which can be used to form an initial guess.

    Hint 2: If you find it difficult to obtain a solution when \(\mu = 100\), then try to increase \(\mu\) in increments from \(0.5\) to \(100\). Each time a new solution is successfully computed, use it as the initial guess for the next (larger) value of \(\mu\). This process of gradually ramping up a parameter and using the previously computed solution as the next initial guess is called natural parameter continuation.

  4. (Optional) Solve the previous exercise using your own implementation of Newton’s method and a built-in algebraic solver (e.g. root or fsolve). Compare the performance of the solver. What solver is more likely to converge?

  5. Extend your code further so that it can solve problems of the form \[ D \tdd{u}{x} + q(x, u, \mu) = 0, \quad u(a) = \gamma_1, \quad \left.\td{u}{x}\right|_{x=b} = \delta. \] As a starting point, set \(q = 1\), \(\gamma_1 = 0\), and \(\delta = 0\). In this case, the exact solution is \[ u(x) = \frac{1}{2D}(x - a)(2b - a - x). \]

  6. (Optional) Try using your code to solve the problem \[ \tdd{u}{x} + 1 = 0, \quad \left.\td{u}{x}\right|_{x=\pm 1} = 0. \] To understand what happens in this case, try solving this problem by hand. Alternatively, integrate the ODE and impose the boundary conditions.

Answers to selected exercises

3. See figure below.

Answer to exercise 3 with \(\mu = 100\).

4. You should see that Newton’s method is more likely to converge, assuming you’re using the default options with the built-in solver.

6. This ODE does not have a solution. Numerical methods will always fail to find a solution.