Coursework 1: MATH 3018/6141 - Numerical methods
Due: Monday 29 November 2021
In this coursework you are to implement the numerical algorithms that are outlined and
described below. The assessment is based on
• correct implementation of the algorithms – 4 marks
• correct use of the algorithms in tests – 3 marks
• figures to illustrate the tests – 2 marks
• documentation of code – 3 marks
• unit testing, robustness and error checking of code – 3 marks.
The deadline is noon, Monday 29 November 2021 (week 9). For late submissions
there is a penalty of 10% of the total marks for the assignment per day after the assign-
ment is due, for up to 5 days. No marks will be obtained for submissions that are later
than 5 days.
Your work must be submitted electronically via Blackboard. Only the Python files
needed to produce the output specified in the tasks below are required.
1
1 Stiff ODEs and transients
The ability to solve an Initial Value Problems (IVP) for y(x), of the form
dy
dx
= f(x,y), y(0) = y0, (1)
depends on the characteristic lengthscales involved in the problem. When the length-
scales differ by many orders of magnitude the system is called stiff and can be very
difficult to solve numerically. In such stiff problems tiny errors in the large scale be-
haviour can swamp the correct transient, small scale behaviour, or vice versa. These
difficulties can usually avoided by using implicit numerical solvers.
An example of such a stiff system is given by the ODE system
d
dx
(
y1
y2
)
=
( −1000y1
1000y1 − y2
)
, y(0) =
(
1
0
)
, (2)
for which the solution is (
y1
y2
)
=
(
e−1000x
1000
999
(
e−x − e−1000x)
)
. (3)
As shown by figure 1, initially y2 jumps very rapidly from 0 to 1 before decaying
very slowly. When such a stiff ODE is solved using a standard explicit algorithm it is
often found that a very small step-length must be take in order to ensure stability of
the result. In other words if too large step length is taken the result of the numerical
method, based on the explicit algorithm, does not even come close to approximating
the solution (because it is unstable). However the very small step length required by
the explicit solver is a waste of computational time and so an implicit solver, which
does not have these stability issues, is more appropriate.
Here we restrict our attention to equations of the form
dy
dx
= Ay + b(x), (4)
where A is an n × n matrix with constant coefficients and the vector b depends only
on the independent variable x (and not on the solution y). Since this equation is linear
its discretisation, using an implicit approach, leads to an algorithm in which a linear
system of equations (i.e. a matrix equation of the form Mu = c) must be solved.
The example above in equation (2), is of the form
A =
(−a1 0
a1 −a2
)
, b = 0. (5)
With the initial data
y0 =
(
1
0
)
(6)
equation (5) has the solution
y =
(
e−a1x
a1
a1−a2 (e
−a2x − e−a1x)
)
. (7)
Stiff behaviour, characterised by an initial rapid transient followed by a slow decay,
occurs for parameter values satisfying a1 a2 > 0.
2
Figure 1: The solution for the system of equation (2) given by
equation (3). Note the different scales required to show the rapid
decay of y1 and the slow decay of y2. Note also the initial rapid
transient behaviour of y2.
2 Algorithms
2.1 Explicit algorithm
You are to implement the standard explict third order Runge-Kutta method, written
RK3. Assume an evenly spaced grid with spacing h (so that xn+1 = xn + h). Then,
given that numerical approximation to y(xn) is yn, the RK3 algorithm used to compute
yn+1 (i.e the numerical approximation to y(xn+1)) from the ODE (4) is
y(1) = yn + h [Ayn + b(xn)] , (8a)
y(2) = 34yn +
1
4y
(1) + 14h
[
Ay(1) + b(xn + h)
]
, (8b)
yn+1 =
1
3yn +
2
3y
(2) + 23h
[
Ay(2) + b(xn + h)
]
. (8c)
where A is the matrix in (4).
2.2 Implicit algorithm
You are also to implement the optimal two stage third order accurate Diagonally Im-
plicit Runge-Kutta method, written DIRK3. Once again assume an evenly spaced grid
with spacing h (so that xn+1 = xn + h). Then, given that numerical approximation to
y(xn) is yn, the DIRK3 algorithm used to compute yn+1 (i.e the numerical approxi-
3
mation to y(xn+1)) from the ODE (4) is
[I − hµA]y(1) = yn + hµb(xn + hµ), (9a)
[I − hµA]y(2) = y(1) + hν
[
Ay(1) + b(xn + hµ)
]
+ hµb(xn + hν + 2hµ), (9b)
yn+1 = (1− λ)yn + λy(2) + hγ
[
Ay(2) + b(xn + hν + 2hµ)
]
, (9c)
where I is the identity matrix, A is the matrix in (4) and the coefficients µ, ν, γ and λ
are defined as
µ = 12
(
1− 1√
3
)
, (10a)
ν = 12
(√
3− 1
)
, (10b)
γ =
3
2
(
3 +
√
3
) , (10c)
λ =
3
(
1 +
√
3
)
2
(
3 +
√
3
) . (10d)
3 Task
1. Implement the explicit RK3 method given in equation (8) as a Python function
x, y = rk3(A, bvector, y0, interval, N)
The input arguments are the matrix A (as defined in (4)), a function bvector,
the n-vector of initial data y0 ≡ y(x0), a list interval giving the start and
end values [x0, xend] on which the solution is to be computed, and N the number
of steps that the algorithm should take (so h = (xend − x0)/N ). The function
bvector should take the form
b = bvector(x)
where the input is the location x (which may vary inside the RK step), and the
output is the vector b as defined in equation (4).
The first output argument x is the locations xj at which the solution is evaluated;
this should be a real vector of length N + 1 covering the required interval. The
second output argument y should be the numerical solution approximated at the
locations xj , which will be an array of size n× (N + 1).
The input to the function should be carefully checked, and the function fully
documented.
2. Implement the implicit DIRK3 algorithm given in equation (9) as a Python func-
tion
x, y = dirk3(A, bvector, y0, interval, N)
4
The input and output arguments follow the same form as for the RK3 algorithm.
Again the function should carefully check its input and should be fully docu-
mented. When solving the linear systems in equations (9a–9b) use the in-built
numpy solvers.
3. Apply your RK3 and DIRK3 algorithms to the system of equations (5) with
the explicit values a1 = 1000, a2 = 1 as in figure 1, using the initial data of
equation (6) over the interval x ∈ [0, 0.1]. This will require defining a function
to provide the (trivial) values of b. Compute the solution for each N where
N = 40k with k = 1, . . . , 10. Compute the 1-norm of the relative error in the
component y2 by comparing to the exact solution in equation (7) (for x 6= 0) as
‖Error‖1 = h
N∑
j=2
∣∣∣∣ (y2)j − ((y2)exact)j((y2)exact)j
∣∣∣∣ . (11)
Plot the error against h, using an appropriate scale. By fitting an appropriate
curve to the data using numpy.polyfit, which should also be plotted, give
evidence to show that the algorithm is converging at third order. In addition plot
the solution computed with the highest resolution and the exact solution against
x on appropriate scales (in one figure for each algorithm, but different variables
in subplots).
4. Define a new system in which the matrix A is given by
A =
 −1 0 0−99 −100 0
−10 098 9 900 −10 000
 (12)
and the vector function b by
b =
 cos (10x)− 10 sin (10x)199 cos (10x)− 10 sin (10x)
208 cos (10x) + 10 000 sin (10x)
 . (13)
With initial data
y0 =
01
0
 (14)
this system has the exact solution
y =
 cos (10x)− e−xcos (10x) + e−x − e−100x
sin (10x) + 2e−x − e−100x − e−10 000x
 . (15)
Apply both the RK3 and DIRK3 algorithms to this system over the interval x ∈
[0, 1]. Compute the solution for each N where N = 200k with k = 4, . . . , 16.
Again plot the error against h, showing evidence of the convergence rate for
DIRK3 only, and plot the solutions for both algorithms computed at the highest
resolution. All components of y should be plotted against the exact solution,
none on logarithmic scales. When computing the error, only the y3 component
should be used.
5
3.1 Summary and assessment criteria
You should submit the Python code electronically as noted above.
You are expected to submit all the Python code needed to produce the required
output. Ideally there should be a top-level Coursework1.py script that, when run,
produces all output. It is possible to complete the coursework by submitting a single
file, but if you wish to use more files that is fine, provided all are submitted.
The script should produce 7 (seven) plots: for each system (the moderately stiff
case in task 3 and the stiff case in task 4), and for each algorithm (RK3 and DIRK3),
two plots are required. The first plot should show the behaviour of the errors against
h and should be annotated to show the convergence rate where appropriate (this plot is
not required for the RK3 algorithm in the stiff case). The second plot should show the
exact and numerical solutions for each component on appropriate scales in individual
subplots (e.g., as in figure 1). All plots should be suitably labelled and clear.
The primary assessment criteria will be the correct implementation of the RK3
and DIRK3 algorithms. The correctness of the algorithms must be clear, and must be
demonstrated by completing the tests shown.
The secondary assessment criteria will be the clarity and robustness of your code.
Your functions and scripts must be well documented with appropriate docstrings and
internal comments describing inputs, outputs, what the function does and how it does
it. The code should also be well structured to make it easy to follow. Input must
be checked and sensible error messages given when problems are encountered. The
clarity of the output (such as plots or results printed to the command window) is also
important.
Code efficiency is not important unless the algorithm takes an exceptional amount
of time to run.
For those interested in applications of these techniques to current problems, the
review of Gottlieb et al. includes a broad range of examples, and the paper of Pareschi
and Russo gives a specific example where DIRK methods are important.
6