Assignment 3

General Feedback

The assignment was solved well in general but here are some common problems that occurred.

  • many people used a linear scale to generate the data for the error plot. If you are going to plot on a log scale it is more efficient to use points equidistant on a log scale. Some people did not make a log-log plot.

  • many people did not update the either the t or N after every step

  • some did not realise that one needs to use a fixed time interval to evolve the ODE over (for example from t=0 to t=1000) and therefore one has to change the size of dt like dt = (t1 - t0) / n_steps to adapt to changing the number of steps. The time used for the analytical reference value should also match!

  • one elegant way of accessing the last element in an array or list a is to use -1 as the index.

a = [ 1, 2, 3, 4, 5, 6, 7, 8]
last_element = a[-1]

  • always use the functions you had written, and importantly tested, at the start. Don’t copy the code further down the notebook. At best this is inefficient, at worst if can lead to mistakes.

  • Some students were still missing labels, title or a legend, or labelling axes as log when a log scale had been used (you should only put log in the label if you have taken and plotted the logs)

Model solution

Automatically marked cells:

# define a function to calculate the mean lifetime from the half life
def meanLifetime(halfLife):
    mean = halfLife/ numpy.log(2)
    return mean

T_HALF = 5730.
TAU = meanLifetime(T_HALF)


def f_rad(N, t):
    return - N / TAU

def analytic(N0, t):
     return N0 * numpy.exp(-t/TAU)

def solve_euler(f, n0, t0, dt, n_panels):
    # Initialise variables
    n = n0
    t = t0
    n_t = numpy.empty((n_panels+1))  # Array ready for updated counts
    
    # For loop generating decay data in n_panels steps of width dt using Euler's method
    for i in range(n_panels+1):
        n_t[i] = n
        n += dt*f(n, t)
        t += i*dt  # For completeness, but this ODE does not need time.
    return n_t

def solve_RK4(f, n0, t0, dt, nsteps):
    res_n = numpy.array([0.0] * ( nsteps + 1 ))

    current_n = n0
    current_t = t0
    
    res_n[0] = current_n
    
    for i in range(1,nsteps+1):
        
        k1 = f(current_n          , current_t)
        n1 = current_n + k1* dt/2.0
        k2 = f(n1, current_t + dt/2.0)
        n2 = current_n + k2* dt/2.0
        k3 = f(n2, current_t + dt/2.0)
        n3 = current_n + k3* dt
        k4 = f(n3, current_t + dt )
        k = 1.0/6.0 * ( k1 + 2 * k2 + 2 * k3 + k4 )
        current_n += dt * k
        current_t += dt
        res_n[i] = current_n
    return res_n

Plotting task

This is an examples for the plotting task.

n0 = 1000
tmax = 3000.
target = analytic(n0, tmax )
err_euler = []
err_rk4 = []
ntest = [ 2**expo for expo in range(20) ] 
for n_step in ntest:    
    res_euler = solve_euler(f_rad,n0, 0, tmax/n_step, n_step)
    res_rk4 = solve_RK4(f_rad,n0, 0, tmax/n_step, n_step)
    err_euler.append( abs( res_euler[-1] - target) )
    err_rk4.append( abs( res_rk4[-1] - target) )
    
plt.plot(ntest, err_euler, label='Euler')    
plt.plot(ntest, err_rk4, label='RK4')    
plt.loglog()
plt.legend()
plt.xlabel('Number of steps $N$')
plt.ylabel('$|n_{num}(t=3000)-n_{ana}(t=3000)|$')

png

Using scipy

Again if we really had to solve this problem we would use scipy. The solve_ivp function provides access to a range of different methods. The default is called RK54 which uses the 5th order Runge-Kutta method to compute the values together with the 4th order to estimate the error. It needs the function for the derivative to have both the value and the time, this is why ours had this, although in the opposite order.

import scipy.integrate as integrate
def f(t,n) :
    return f_rad(n,t)

output = integrate.solve_ivp(f,(0,tmax),numpy.array([n0]),rtol=1e-9)
print((output.y[0][-1]-target)/target)

1.3500933701583496e-10