Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
MATH 3018/6141 - Numerical methods
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).