Navigation :

Writing More Complex Tests

Convergence

Because of the ubiquity of differential equations in physical problems, it’s a common exercise to calculate the numerical derivatives of functions via a Central difference formula:

$$ f’(x) = \frac{f(x + h) - f(x - h)}{2h} + \mathcal{O}(h^2) $$

How can we test this? We’ll decide to call our function CentralDiff. First we can think about simple tests that check that an error is raised if input arguments do not make sense. We know what parameters our function will take, so we define an ‘interface’ to the function:

def CentralDiff(func, x, h):
    """
    Code to be implemented
    """
    pass

We’ll need to pass a function, the x position where the derivative is computed and the step size. Tests of the input arguments are:

import numpy as np
import pytest

def test_inputs_error_raised():
    # It does not make sense for h to be exactly zero, so this should raise
    # an error.
    with pytest.raises(ValueError):
        CentralDiff(np.sin, 0.0, 0.0)


    # If an argument other than a function is passed, a TypeEror should be raised.
    with pytest.raises(TypeError):
        CentralDiff(2.0, 0.0, 1e-5)

    # We check that non-numeric inputs raise a TypeError
    with pytest.raises(TypeError):
        CentralDiff(np.sin, '0.0', 1e-5)

    with pytest.raises(TypeError):
        CentralDiff(np.sin, 0.0, '1e-5')

Save the non-implemented function to a file called ‘functions.py’, and save the test to test_functions.py

We can test the function against known derivatives of analytic functions that we calculate by hand.

Exercise: Convergence

How would you test the convergence? Write a test which checks the error on the derivative of Sin at a known point for different values of h. What can you check numerically? Write an implementation of the numerical function.

If you struggle, it may help first to implement the CentralDiff function and look at the results, and plot the error.

Exercise: Newton Raphson

The Newton-Raphson algorithm finds minima of functions with the following recursion relation:

$$ x{n+1} = x{n} - \frac{f(x_n)}{f’(x_n)} $$

We normally check convergence by looking at the difference between the last step and the previous step, and stopping when we reach a tolerance for $$ abs(f(x_n)) < \epsilon $$

Implement the algorithm as a function. Write a test that checks this works. {% /panel %}

Exercise: Euler ODE Solver

How would one go about testing an ODE solver? Consider the prototypical simple 2nd order ODE - the simple harmonic oscillator:

$$ m \frac{\mathrm{d}^2 x}{dt^2} + k x = 0 $$

We know from A Level maths that the general solution to this is: $$ x(t) = A \cos{\omega t} + B \sin{\omega t} $$ with $$\omega = \sqrt{k / m} $$

To compute the solution numerically, we have to decompose the ODE:

$$\frac{\mathrm{d}v}{\mathrm{d}t} = \frac{-kx}{m} $$ $$\frac{\mathrm{d}x}{\mathrm{d}t} = v $$

Write a very simple Euler integrator for this problem. Use the initial conditions:

$$ x(t = 0) = 1.0 $$ $$ v(t = 0) = 0.0 $$

with paramters (k = 1) and (m = 1)

Think about how you would design the interface to this so that it can be generally applied to different problems. We know the analytic solution for this problem, so play around with how the error on the position of a mass increases as you uses. You need to have at least two functions to unit test here!

Hint: Think about both the solution and properties you can compute about the system you’re simulating!