EMAT20920: Numerical Methods in MATLAB
Question 1: Root-finding
(a) For each of the equations (i)–(iii) below, state
• the number of solutions of the equation in the given range,
• a bracket for each of the solutions,
• the values of the solutions, as accurately as possible (using a built-in MATLAB root-finding
You should report not just the solutions, but also the code used to find them, and any supporting
material necessary to justify your solutions.
x−1 = tanh(2x − 1) for x ∈ R,
for −2π < x < 2π,
5 + 1000x
4 + 35460x
3 + 286290x
2 − 332910x − 39340
4 − 299
= 499 for x ∈ R,
(b) The bisection method — which finds a root of a continuous function f(x) in an interval [a, b] on
which f changes sign — can be modified by using a point other than the midpoint m of the interval
as the estimate of the root. One alternative is to instead use the point c where the straight line
joining the points (a, f(a)) and (b, f(b)) crosses the x-axis; that is
af(b) − bf(a)
f(b) − f(a)
EMAT20920 2019-20 1 Coursework Assessment
Use this modified bisection method (which uses c in equation (1), instead of the midpoint) to find
the root of f(x) = cos(x) − x in the interval [0, 1], accurate to 13 decimal places. Use the results
to plot a graph that shows the order of convergence of the method, and use the graph to find the
(c) The iteration
xn−2f(xn−1) − xn−1f(xn−2)
f(xn−1) − f(xn−2)
can be used, as n → ∞, to find a root x = limn→∞
xn of the function f, given initial values x0 and x1.
Use equation (2) to find the root of f(x) = cos(x) − x, given initial values x0 = 0 and x1 = 1,
accurate to 13 decimal places. Use the results to plot a graph that shows the order of convergence
of the method, and use the graph to find the order.
Question 2: Numerical calculus
(a) The expression
b − a
a + b
b − a
a + b
a + b
b − a
can be used to approximate the integral R b
Use this definition to create a MATLAB function that calculates a numerical approximation to
an integral R β
f(x) dx — given an integrand f, region [α, β], and absolute error tolerance —
using a composite version of equation (3); i.e., where the interval [α, β] is divided into equal-width
subintervals, on each of which we use equation (3). Use your function to find the degree of accuracy
and order of the composite rule, explaining your reasons in each case.
(b) The continuous random variable X has probability density function
log(1 + kx) e−x/k x > 0
where k > 0 is a constant, and C(k) = 1/
log(1 + kx) e−x/k dx. Use MATLAB to plot a
(single) graph of the mean, median, and mode of X, as a function of k.
EMAT20920 2019-20 2 Coursework Assessment
(c) The expressions
−f(a + 2h) + 8f(a + h) − 8f(a − h) + f(a − 2h)
−11f(a) + 18f(a + h) − 9f(a + 2h) + 2f(a + 3h)
can both be used, as h → 0, to approximate the derivative of a function f(x) at x = a.
Use equations (4) and (5) to approximate the derivative of f(x) = x sin(1/x) at x = 3/π, and plot
a graph of the error between the approximate and true value of the derivative in both cases, for a
range of values of h sufficient to illustrate the different sources of error in the formulae. Describe
the results that you find, and explain which is the better formula and why.
Question 3: Ordinary differential equations
The following set of ODEs is a simple model of infectious disease dynamics in a population:
= −βSI, (6)
= βSI − µI, (7)
= µI. (8)
The variables represent the proportions of the population at time t who are 1. healthy but susceptible
to the disease, S(t), 2. infected with the disease, I(t), and 3. recovered from the disease, R(t). Recovery
from the disease confers immunity, and there is no death (either from the disease, or from natural causes).
The parameter β is the infection rate, and µ is the recovery rate.
(a) Mathematically (without using MATLAB) show that the total population S + I + R is a conserved
quantity, i.e. that it remains fixed for all time. Explain why that means we can set S + I + R = 1,
and why we only need to integrate equations (6) and (7) to solve the SIR system.
(b) Suppose that at the intial time t = 0, 1% of the population is infected with the disease and everyone
else is susceptible, so S(0) = 0.99 and I(0) = 0.01. Set β = 2 and µ = 1. Use the forward Euler
method with stepsize h = 0.1 to solve the SIR model up to time t = 15. Plot a graph of S(t), I(t),
R(t) against t (in a single graph), with nice-looking labels, distinguishable lines, etc., and explain
what you see.
EMAT20920 2019-20 3 Coursework Assessment
(c) Write MATLAB code that uses both the forward Euler and the fourth-order Runge-Kutta methods
to solve the SIR model, with the same parameter values and initial conditions as in part (b) above,
for a given stepsize h. Use it to plot a single graph of the final R(t) value against stepsize h, for
different values of h between 100 and 10−3
. Explain what you see.
(d) Use an in-built MATLAB ODE solver to solve the SIR model, with the same parameter values and
initial conditions as in part (b) above, and find the time for which I(t) is a maximum.
(e) Fix the recovery rate µ = 1. Solve the SIR model, using an in-built MATLAB ODE solver, until
the dynamics reaches a steady state, for various values of β, and plot a graph of the final R(t) value
against β. Explain the graph you have obtained.