SEMT20001-R: Principles of Computational Modelling
COURSEWORK REASSESSMENT: QUESTION 4
Question 4: Equation-based modelling
(a) (i) Write Python code that numerically solves an ODE system
(1)
given a right-hand-side function f, timespan [t0, tf], initial condition y0, and number of steps N / step-size h = (tf - t0)/N, using the time-stepping algorithm
(2)
where
(2 marks)
(ii) Using both the forward Euler method, and the method defined in part (a)(i), solve the ODE initial value problem
(3)
for t ∈ [1, 3], with step-size h = 0.5. Plot a graph of the numerical solutions, together with the exact solution
(4)
and comment on the results. (4 marks)
(iii) Plot a graph of the absolute error of the numerical solution of equation (3), using the method given in part (a)(i), versus step-size h, for a suitably-chosen range of h. You should use the Euclidean distance between the numerical solution
(5)
to define the error. Use the graph to find the order of the truncation error as h → 0. (4 marks)
(b) The ODE system dt/dy = f(t, y), written in terms of components y = (S, E, I, R),
represents an equation-based model of the transmission of a non-fatal disease with four states:
in which the population is divided into four groups:
• Susceptible: healthy and not yet infected with the disease,
• Exposed: infected with the disease, but not yet infectious,
• Infectious: both infected and infectious (i.e., able to pass on the disease to others),
• Recovered: healthy, with long-lasting immunity.
The variables in the model, S(t), E(t), I(t), and R(t), are the (non-dimensional) proportions of the total population in the four states at any time t. The parameters are b the birth rate, d the natural-causes death rate (assumed equal across the population groups), β the infection rate, γ the recovery rate, and α the infected to infectious rate (so 1/α is the average latency period). All the parameters are assumed to be constant and positive, and measured per simulated day.
(i) Use solve_ivp to solve the SEIR model (6), with parameter values b = d = 0.01, α = 0.25, β = 0.2 and γ = 0.1, for t ∈ [0, 3 × 365], with initial conditions S(0) = 0.99, E(0) = 0.01, I(0) = 0 and R(0) = 0, and plot a graph of the solution versus time. Explain what the results mean in the context of the infection in the population.
How long does it take for the solution to reach a steady state, where the magnitude of the rate of change of the solution, |dt/dy|, is less than 10-5 ? (4 marks)
(ii) Plot a graph of the steady-state solution of the SEIR model (6) versus β, for β ∈ [0, 0.5], by using solve_ivp to solve the SEIR model (6) until the dynamics reach a steady state for a range of values of β . The other parameters should be fixed to the values given in part (b)(i).
Explain the graph you have obtained, in the context of the infection in the population. (4 marks)
(c) As described in the lecture notes, the Robertson system is an equation-based model of a reaction between three chemicals X , Y , Z, with concentrations x(t), y(t), z(t):
where the reaction rates are k1 = 0.04, k2 = 3 × 107 , k3 = 104 .
Write a Python function that measures the performance of a solve_ivp solver method, when solving the Robertson system (7) for t ∈ [0, 5] with initial conditions x(0) = 1, y(0) = z(0) = 0, by returning:
• the number of right-hand-side function evaluations of solve_ivp,
• the average run time of solve_ivp,
given a solve_ivp method (the lectures used method= ' RK45 ' and method= ' Radau ' ), and desired relative and absolute tolerances, rtol and atol.
Use your code to plot graphs of
• the number of right-hand-side function evaluations of solve_ivp,
• the average run time of solve_ivp,
versus relative error tolerance rtol, for at least four solve_ivp solver methods, with values of rtol between 10-4 and 10-8 and atol=1e-3*rtol.
Interpret the results, and explain how they might influence your choice of method when using solve_ivp. (7 marks)