Question 1:    Errors and orders

(a)  The iteration

xn =                                                                                               (1)

noo

Write MATLAB code that uses equation (1) to find the root of f (x) = e_x -x, accurate to 13 decimal places, given initial values x0  = 0 and x1  = 1. Use the results to plot a graph that shows the order of convergence of the method, and use the graph to find the order.

(6 marks)

(b)  The expressions

-f (a) + 3f (a + h) - 3f (a + 2h) + f (a + 3h) h3                                                 ,

and

-f (a - 2h) + 2f (a - h) - 2f (a + h) + f (a + 2h) 2h3

can both be used, as h → 0, to approximate the third derivative of a function f (x) at x = a.

Write MATLAB code that uses equations (2) and (3) to approximate the third derivative of f (x) = x cos(1/x) at x = -2/π, and use it to plot a graph of the absolute error between the approximate and true value of the derivative in both cases.  The range of values of h should be sufficient to illustrate all the different sources of error in the formulae.

What are the orders of the errors, in each case? Which is the better formula, and why?

(7 marks)

(c)  The iteration

yk+1 = yk + hf ╱tk +  , yk + f (tk, yk)、

tk+1 = tk + h

can be used to find a numerical solution {(tk, yk), k = 0, . . . , N } of the first-order ODE y˙ = f (t, y), given a timespan t e [t0 , tf], initial condition y0  = y(t0 ), and stepsize h = (tf - t0 )/N for some positive integer N .

Write MATLAB code that uses equation (4) to find the solution of

dy

and use it to plot a graph of the global truncation error as a function of h. Use the graph to find the order of the global truncation error.

What other numerical errors might we expect to see when using (4) to solve an ODE?

(7 marks)

Question 2:    Numerical integration

(a) Use built-in MATLAB commands to find the values of the following integrals (you should report

both the command, and the value of the integral):

1       cos x                                                           o                   1

_1  ^1 + x2                                                                 _o  (1 + IxI) 尸IxI

(2 marks)

(b)  The continuous random variable X has probability density function (PDF)

fX (北) = ,．(0(α) eα+x_α ez       o(北)

where α > 0 is a constant. Use MATLAB to plot four graphs, as follows:

(i) the PDF as a function of 北, for 0 < 北 < 3 and α e {0.5, 1, 1.5, 2},

(ii) the cumulative distribution function  (CDF) as a function of 北,  for 0  <  北  <  3 and α  e

{0.5, 1, 1.5, 2},

(iii) the mean and mode of X, as functions of α, for 0.1 < α < 2,

(iv) the standard deviation and interquartile range of X, as functions of α, for 0.1 < α < 2. In each case, you should show the commands you use, as well as the final graphs.

(8 marks)

(c) Write a numerical integration function in MATLAB, that will find a numerical approximation (without using any of the built-in numerical integration commands) to the integral

b

f (北)d北

a

accurate to within a user-specified tolerance. As a minimum, your function should have as inputs the function f to be integrated, the domain of integration  [a, b], and the desired tolerance, and should return the approximation to the integral found using one of the composite Newton-Cotes rules.  It should also be able to deal with the case that the approximation does not converge (for example, if the integral is unbounded).

Other possible features that you might choose to add to your code could include:

● allowing choice of relative and/or absolute error tolerances,

● display of diagnostic information during or after the progress of the algorithm,

● help text explaining usage,

● checks on the input parameters,

● implementation of a higher-order method (e.g. composite Gaussian quadrature).

Together with the function, you should write a brief description of how your code works, and details of any tests you conducted to make sure it works correctly.

(10 marks)

Question 3:    Numerical solution of ODEs

The dynamics of a rock-scissors-paper type game, played amongst a large group of people, are modelled by the ODE system

= αxy - γxz,

dt

dy

= γzx - βzy,

dt

where x(t), y(t) and z(t) are the proportions of the population at time t who are playing rock, scissors and paper, respectively, and α > 0, β > 0, γ > 0 are constants.

The αxy and -αyx terms on the right-hand side of equations (5) and (6), respectively, represent the fact that rock (x) beats scissors (y): when they meet (which occurs with frequency αxy per unit time), the

proportion of the population playing rock increases at the same rate as that playing scissors decreases. (a) Explain how the other types interact by referring to the ODEs, equations (5)– (7).

(1 mark)

(b)  Show algebraically  (i.e., without using MATLAB) that (i) x + y + z, and (ii) xyz are conserved

quantities  (i.e., that they do not vary with time) if α = β = γ .  Which, if either, is a conserved quantity if α , β , γ are not all equal?

(3 marks)

(c) Write MATLAB code that uses the forward Euler method to solve the  ODE system given by equations  (5)– (7).   Use your code to solve the system with  α  =  β  = γ  =  1,  initial condition (x, y, z) = (0.5, 0.3, 0.2) at t = 0, from t = 0 to t = 200, with

(i) timestep h = 0.02, and (ii) timestep h = 0.002.

For each of (i) and (ii), plot graphs of x, y , z , x + y + z, and xyz against time t.

Explain what you see in words. Are the results consistent with the mathematical ones obtained in part (b)? Why/why not? What is the effect of changing h?

(7 marks)

(d) Write MATLAB code to solve the ODEs (5)– (7) using 。de45 , with the same parameter values, initial conditions and total integration time as in part (c), making sure that the conserved quantities remain constant. Plot the results (x, y , z , x + y + z, and xyz against time t), as in part (c).

(3 marks)

(e) Write MATLAB code that uses 。de45 to find both the peak and trough amplitude of the periodic

solutions x, y, and z, at the end of the simulation, as shown for x in the figure below.

160                170                180                190                200

t

Use your code to explore what happens to the solutions when the parameter α varies between 0.25  and 4,  by plotting the peak  and trough  amplitudes of x,  y,  and  z,  versus  α .   Keep the other parameters, initial conditions, and total integration time fixed (i.e., β = V = 1, (x, y, z) = (0.5, 0.3, 0.2) at t = 0, solved from t = 0 to t = 200).

(6 marks)