Solving Ordinary Differential Equations (ODEs)

You can download the slides for this lecture here

The intentions behind this exercise are:

To give you a clear goal to meet
To give you a chance to revise or learn a few things. Specifically:
- The abstract base class pattern;
- Unit testing with unittest;
- Ordinary differential equation solver implementations.

A nice feature of this assignment is that you are provided with a testing infrastructure which gives an implicit specification of the code you are required to write. This means that when you've written code such that all the tests pass then you will know that you have completed the exercise.

Getting the code for these exercises

Download the file ODEsAM.zip. You have been provided with one source class (AbstractOdeSolver) and with two test files. You will not need to edit any of these 3 source files. If you do make edits to them during the course of the exercise then make sure that your final version works with the original versions of these files.

Forward Euler

Running first test:

Check to see what running the simplest test-suite TestOdeSolvers.py does at this stage. It should stop with ModuleNotFoundError: No module named 'ForwardEulerOdeSolver'. In order for the test to run, you need to create the ForwardEulerOdeSolver class in ForwardEulerOdeSolver.py. This class should inherit from AbstractOdeSolver, and should provide the one method: Solve(). (This method's implementation does not yet have to contain any code other than pass.) Now check that the TestOdeSolvers.py test-suite runs. It is expected to throw exceptions at this stage.

Writing a forward Euler solver:

The good news is that 1 out of 3 of the tests in the TestOdeSolvers test-suite pass immediately. The first test test_abstract_class_methods verifies that all the functionality I have written in the abstract base class works and consequently the test passes. The other two tests which check ForwardEulerOdeSolver.Solve() works correctly are the tests which are failing. If you need a revision tutorial about the Forward Euler method then please ask now.

Your task is to implement the Solve() method so that it solves from the initial values over a time-interval with a fixed time step using forward Euler. On exit from the method you should have populated solutionTrace with a list of

(x,y)

values (as pairs/tuples) and timeTrace with the list of corresponding times. The test-suite uses two simple model ordinary differential 2-d systems of equations. These are plugged into the solver as right-hand side functions. They are defined at the top of the test-suite file and both have straightforward analytic solutions.

As you develop your Solve() method you should reduce the number of failing lines. One or two of the test lines are harder to satisfy than others. Attempt to have no failures but the fewer the better. Please make sure that you have comments in your code.

Here's a version which first builds the lists to the correct length and then writes in using indexing.

class ForwardEulerOdeSolver(AbstractOdeSolver):
    def Solve(self):
        """ The forward Euler solve method """
        # Fulfil the "throw if not setup" test
        if self.numberOfTimeSteps <= 0 :
            raise Exception("Timestepper not set up")

        # There's a test which runs the Solve method twice
        # and checks that the data from the first run is erased
        self.solutionTrace=[None]*(self.numberOfTimeSteps+1)
        self.timeTrace=[0]*(self.numberOfTimeSteps+1)

        # Initial state is recorded
        self.solutionTrace[0] = self.initialValues
        self.timeTrace[0] = self.startTime

        # Set up loop variables
        time = self.startTime
        values = self.initialValues

        for i in range(0, self.numberOfTimeSteps):
            # Euler step
            dvdt = self.rhsFunction(values, time)
            x = values[0] + self.timeStepSize*dvdt[0]
            y = values[1] + self.timeStepSize*dvdt[1]

            # Update solution and time to new values
            values = (x,y)
            time = self.startTime + (i+1)*self.timeStepSize
            self.solutionTrace[i+1] = values
            self.timeTrace[i+1] = time
        # Initial state is recorded
        self.solutionTrace.append(self.initialValues)
        self.timeTrace.append(self.startTime)
        # Set up loop variables
        time = self.startTime
        accumulated_time = self.startTime  # Wrong way to do this
        values = self.initialValues

        for i in range(0, self.numberOfTimeSteps):
            # Update solution and time to new values
            values = (x,y)
            time = self.startTime + (i+1)*self.timeStepSize
            accumulated_time += self.timeStepSize # Wrong way to do this
            self.solutionTrace.append(values)
            self.timeTrace.append(accumulated_time)

Here's an outline of a version with appends to the back of the lists. This code also illustrates where testing may fail if the time-stepper uses repeated addition.

A Higher-Order Solver

Now run TestHigherOrderOdeSolverRunner. This new test-suite requires you make another class HigherOrderOdeSolver.py which should again provide only a Solve() method. It tests the convergence behaviour of your two solvers and checks that the higher-order one is better. My model answer (which passes the test) uses a standard second order Runge-Kutta method. Again work on your solution until there are few lines (ideally no lines) of failing tests.

class HigherOrderOdeSolver(AbstractOdeSolver):
    def Solve(self):
        """This is a Runge-Kutta (RK2) solve method"""
        # Fulfil the "throw if not setup" test
        if self.numberOfTimeSteps <= 0 :
            raise Exception("Timestepper not set up")

        # There's a test which runs the Solve method twice
        # and checks that the data from the first run is erased
        self.solutionTrace=[None]*(self.numberOfTimeSteps+1)
        self.timeTrace=[0]*(self.numberOfTimeSteps+1)

        # Initial state is recorded
        self.solutionTrace[0] = self.initialValues
        self.timeTrace[0] = self.startTime

        # Set up loop variables
        time = self.startTime
        values = self.initialValues

        for i in range(0, self.numberOfTimeSteps):
            #  Use gradient to get halfway point
            dvdt = self.rhsFunction(values, time)
            halfwayx = values[0] + self.timeStepSize*dvdt[0]/2.0
            halfwayy = values[1] + self.timeStepSize*dvdt[1]/2.0

            # Recompute gradient at halfway point
            dvdt = self.rhsFunction((halfwayx, halfwayy),
                                    time+(self.timeStepSize/2.0))

            # Use this new gradient to update from the original value
            x = values[0] + self.timeStepSize*dvdt[0]
            y = values[1] + self.timeStepSize*dvdt[1]

            # Update solution and time to new values
            values = (x,y)
            time = self.startTime + (i+1)*self.timeStepSize
            self.solutionTrace[i+1] = values
            self.timeTrace[i+1] = time

Using a library function

Write a solver implementation which uses scipy.integrate.odeint at its core.
You should be able to do this with very few lines of code, including just a single call to odeint. Test it in a similar way to the HigherOrderOdeSolver, using progressively smaller time-steps.

What do you notice about the reported error?
What control do you have over the error?

import numpy as np
from scipy.integrate import odeint
class OdeIntOdeSolver(AbstractOdeSolver):
    def Solve(self):
        """ The solve method using scipy.integrate.odeint"""
        # Fulfil the "throw if not setup" test
        if self.numberOfTimeSteps <= 0 :
            raise Exception("Timestepper not set up")
        stopT = self.startTime+self.numberOfTimeSteps*self.timeStepSize
        times = np.linspace(self.startTime,stopT,self.numberOfTimeSteps+1)
        self.timeTrace = list(times)
        self.solutionTrace =\
            odeint(self.rhsFunction, self.initialValues, times)

With the default setting and the same test as before I get something similar to:

#step_size         euler_error        hi_error    oi_error
6.283185307179586  4.442882938158366  14.647777676841754 1.0706378471935034e-07
3.141592653589793  7.088890675779233  14.73085508342019  1.1477237333346198e-07
1.5707963267948966 6.5841823894717395 2.953368442233411  1.1211067003602387e-07
0.7853981633974483 3.103938595129146  0.4282090775868214 1.1392372439136046e-07

This is because odeint has absolute and relative error tolerances set and it uses adaptive time-steps between the time-steps which we care about. That means that the first

2\pi

step around the circle is about as accurate as the finer tests which follow it.

We can control the error by setting the error tolerances, with

 self.solutionTrace =\
     odeint(self.rhsFunction, self.initialValues, times, rtol=1e-12, atol=1e-12)

we get:

#step_size         euler_error        hi_error    oi_error
6.283185307179586  4.442882938158366  14.647777676841754 2.4189691882833333e-12
3.141592653589793  7.088890675779233  14.73085508342019  2.2407390802495764e-12
1.5707963267948966 6.5841823894717395 2.953368442233411  1.990068084042871e-12
0.7853981633974483 3.103938595129146  0.428209077586821  1.6920994988211508e-12

Convergence Behaviour

[Extension] Use the test output in TestHigherOrderOdeSolver as the model for plotting convergence behaviour for these two solvers and at least one other solver that you should write. Use matplotlib or some other plotting program to make a log-log plot with time-step on the x-axis and error on the y-axis.

C++

[Extension] If you love C++ then you are welcome to investigate coding ODE solvers in C++ by investigating the assignment which formed the basis of this exercise: https://www.cs.ox.ac.uk/people/joe.pitt-francis/oxfordonly/FridayAssignment.pdf

ODE Solvers - AM Exercises