Department of Engineering Mathematics
EMAT20920: Numerical Methods in MATLAB
WORKSHEET 5
COURSEWORK ASSESSMENT
The following questions together form the assessment for the unit.
There are 3 questions, each of which carries 20 marks. Answer all questions.
The maximum mark for this assessment is 60 marks.
Your answers should be submitted on Blackboard by 23:59 on 13 December 2019.
You must submit both (1) a single self-contained document (.pdf or .docx), that reports not only
your answers to the questions, but also explains how you used MATLAB to find them; M-files,
where necessary, should be included in this document (using appendices, if the code is extensive),
and (2) a single .zip file containing all M-files.
Your submission must be your own individual effort. Only answers contained in (1), and M-files
contained in both (1) and (2) will be marked. Your submission must clearly state your student
number.
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
routine).
You should report not just the solutions, but also the code used to find them, and any supporting
material necessary to justify your solutions.
(i) x2 ex−1 = tanh(2x− 1) for x ∈ R,
(ii) tan
(
25/2x
3
)
= tan
(
2x
3
)
for −2pi < x < 2pi,
(iii)
50000x5 + 1000x4 + 35460x3 + 286290x2 − 332910x− 39340
300x4 − 299 = 499 for x ∈ R,
(10 marks)
(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
c =
af(b)− bf(a)
f(b)− f(a) . (1)
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
order.
(5 marks)
(c) The iteration
xn =
xn−2f(xn−1)− xn−1f(xn−2)
f(xn−1)− f(xn−2) (2)
can be used, as n→∞, to find a root x = lim
n→∞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.
(5 marks)
Question 2: Numerical calculus
(a) The expression
b− a
18
[
5f
(
a+ b
2
−
√
3
5
b− a
2
)
+ 8f
(
a+ b
2
)
+ 5f
(
a+ b
2
+
√
3
5
b− a
2
)]
(3)
can be used to approximate the integral
∫ b
a f(x) dx.
Use this definition to create a MATLAB function that calculates a numerical approximation to
an integral
∫ β
α 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.
(7 marks)
(b) The continuous random variable X has probability density function
fX(x) =
C(k)x2 log(1 + kx) e−x/k x > 00 otherwise
where k > 0 is a constant, and C(k) = 1/
∫∞
0 x
2 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.
(6 marks)
EMAT20920 2019-20 2 Coursework Assessment
(c) The expressions
−f(a+ 2h) + 8f(a+ h)− 8f(a− h) + f(a− 2h)
12h
, (4)
and
−11f(a) + 18f(a+ h)− 9f(a+ 2h) + 2f(a+ 3h)
6h
(5)
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/pi, 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.
(7 marks)
Question 3: Ordinary differential equations
The following set of ODEs is a simple model of infectious disease dynamics in a population:
dS
dt
= −βSI, (6)
dI
dt
= βSI − µI, (7)
dR
dt
= µ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.
(2 marks)
(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.
(5 marks)
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.
(5 marks)
(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.
(3 marks)
(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.
(5 marks)
EMAT20920 2019-20 4 Coursework Assessment