Question 1 (From Chapter 1 of the Numerical Analysis part of the Text)

The diagram above shows two corridors of widths a and b (a ≤ b) which meet at right angles. The longest rectangular object (desk or piano etc) of width w (w < a) that can be manoeuvred around the corner has a length l, given by:

l = (mb + a -w√(1+m²))√(1+1/m²)                         (1)

provided that the object cannot be moved vertically. The value of m can be found from the condition that in the interval [0, (a/b)^1/3 ] m is the only root of the equation:

(b² - w² )m⁶ + w²m⁴ - 2abm³ + w²m² + a² - w² = 0      (2)

(a) Write a program using MATLAB to solve Eq. (2) by the Regula Falsi method, and hence find m and l. As a test, for a = 50, b = 70, and w = 40, l = 87.4723. Discuss how the initial points should be chosen.

(b) For thin objects like a ladder, w ≈ 0. When w = 0 show that the method will not work if the interval [0, 1] is searched. Explain why this happens. Hint: solve Eq. (2) analytically for w = 0.

Question 2 : a) (from Chapter 2 of the Numerical Analysis part of the Text) Write a program that implements the LU factorisation algorithm 2.3, with partial pivoting, given in the Text and in the supplementary section with MATLAB versions of the algorithms. Your code should have separate functions for the factorization, forward and back-substitution steps. Check your code by analysing the problem solved on page 41 of the Text:

b) (From Chapter 3 of the Numerical Analysis part of the Text) Implement the least squares fitting algorithm 3.2 given in the Lecture Notes and in the supplementary section with MATLAB versions of the algorithms. Use your LU factorisation program (from part (a) of this question) to solve the linear equations and a separate function to evaluate the various expressions. Check your program by analysing the problem given on pages 67-68 of the Text.

c) As an example for an application, you are given the following x-y data representing student attendance at lectures in a particular course:

Use your program to find the coefficients for a least squares fit to these data using the expansion:

f(x) = a₀ + a₁x + a₂x² + a₃x³

Write your program in MATLAB so that it prints the value of f (x) at each of the x data points. How well does f(x) fit the data?

Question 3:

The SIR model - where "S" stands for susceptible, "I" for infected, and "R" for recovered - is one of the simplest mathematical models for the spread of an infectious disease among a population. The equations are:

dS/dt = -aS(t)I(t)

dI/dt = aS(t)I(t) - bI(t)

and

R(t) = N_tot - S (t) - I (t)

where t is time, N_tot is the total number of people in the population, a and b are constants. The constant b is related to the period of infection, p, by p = 1/b. Thus b = 1/3 means that the average period of infection is 3 days. The model also assumes that recovered individuals do not get re-infected.

Assume that a single infected person with a disease for which a = 0.0001 and b = 1/3 arrives on an island with a population of 10,000, making the total number of people Ntot = 10,001.

(a) Perform one step (with Δt = 1 day) of the numerical integration of these ODE's using the fourth order Runge-Kutta method (5.2) to calculate S (1), I(1) and R(1).

(b) By using the ode45 estimate the number of days required for the disease to ‘run its course', i.e. for I to be < 0.5. (ode45 - intrinsic MATLAB function to solve systems of ODEs.)

Question 4:

The torsion (twisting) of a long prismatic bar is governed by the following Poisson equation:

 ∂²Φ/∂x² + ∂²Φ/∂y² = - 2       (1)

where Φ is the torsion function whose value is zero at the edges of the bar.

To solve (1) for a bar with a rectangular cross-section, of any lengths a in the y-direction and b in the x-direction, the cross-section is to be divided into an equal number (n) internal nodes in both the x- and y- directions. For node (i,j) with i denoting the x-direction, the finite difference approximations to the two derivatives in (1) are:

∂²Φ_i,j/∂x² ≈ Φ_i+1,j - 2Φ_i,j + Φ_i-1,j/h_x²

and

∂²Φ_{i, j}/∂y² ≈ Φ_i,j+1 - 2Φ_i,j + Φ_{i, j-1} x/h²_y           (2)

where h_x is the x-distance between adjacent nodes and hy is the y-distance.

(a) Show that the use of (2) leads to the following equation for Φ_i,j :

Φ_i,j = (Φ_i+1,j + Φ_i-1,j)a'/2 + (Φ_i,j+1 + Φ_i,j-1)b'/2 + a'b'(a² + b²)/(n+1)²      (3)

where a' = a²/a² + b²  and b²/a² + b²

(b) Equation (3) leads to the following form suitable for solution by an iterative technique:

Φ^(k+1)_i,j = Φ^(k)_i,j + resid                                                (4)

where the superscript indicates iteration number. The residual term, resid, is given by

resid = (Φ^(k)_i+1,j + Φ^(k)_i-1,j)a'/2 + (Φ^(k)_i,j+1 + Φ^(k)_i,j-1)b'/2 - Φ^(k)_i,j + a'b'(a² + b²)/(n + 1)²Φ^(k)      (5)

Write a MATLAB program to solve (4) and (5) for user inputs of n, lengths a and b, and maximum number of iterations. The program is to output the maximum value of Φ. As a test run your program for n = 50, a = 62cm and b = 75.5cm. Restrict maximum number of iterations to 1000 and set convergence tolerance, ε = 0.2.

Algorithm for Solving Equation (4) in a rectangular area b x a

In: Number of nodes in both directions, n
Lengths a and b
Maximum number of iterations, max_iter
Convergence tolerance, ε

Out: phi(n, n) - the approximate steady state torsion function

comment: initialise all elements of phi to zero

adash = a^2/(a^2+b^2) bdash = b^2/(a^2+b^2)

last = adashxbdashx(a^2+b^2)/(n+1)^2 do for k = 1, max_iter

maxresid = 0.0 do for i = 2, n+1

do for j = 2, n+1

resid = 0.5x(phi(i+1,j) + phi(i-1,j))xadash + 0.5x(phi(i,j+1) + phi(i,j-1))xbdash - phi(i,j) + last

if |resid| > maxresid then

maxresid = |resid|

endif

phi(i,j) = phi(i,j) + resid

enddo enddo

if |maxresid| < ε then

exit with solution in phi

endif enddo

error: maximum number of iterations exceeded