1 MTHS2008
The University of Nottingham
SCHOOL OF MATHEMATICAL SCIENCES
SPRING SEMESTER 2021-2022
ASSESSED COURSEWORK 2
Submission deadline: Thursday 7th April, 3pm
This coursework contributes 20% towards the overall grade for the module.
This is version 2 of the question paper (14:21 Thursday 17th March, 2022).
Questions are now complete.
Rules:
Each student is to submit their own coursework.
You are allowed to work together and discuss in small groups (2 to 4 people), but you must
write all code by yourself.
You must submit one .py file for each question. A full submissions contains three Python
files: multiple_integration.py, time_steppers.py and finite_difference_1d.py. You are strongly
encouraged to use the Spyder IDE (integrated development environment). Hence you should
not write IPython Notebooks (.ipynb), and you should not use Jupyter.
You may adapt any code that we developed in Python workshops (i.e. you do not have to start
from scratch).
You are allowed and encouraged to use Piazza to obtain clarification of the coursework questions
if needed, as well as ask general Python queries. However, when using Piazza, please ensure
that you do not reveal any answers to others, or directly reveal any code that you have written.
Doing so is considered Academic Misconduct.
Marks breakdown:
This coursework is out of 100 marks:
• Q1 – 30 marks.
• Q2 – 30 marks.
• Q3 – 30 marks.
• Commenting, structure and docstrings– 10 marks.
MTHS2008 Turn Over
2 MTHS2008
Guidelines for submitting the coursework correctly:
For full output data marks, your functions need to return the correct outputs for the particular
input arguments.
For full plot marks, your plot(s) should have the correct data plotted, axis labels, a legend if
appropriate, and a descriptive title. Note that you will not be penalised for not using TeX/LaTeX
fonts.
For full commenting marks, your modules should be commented throughout (where appropriate)
with helpful and suitably succinct comments.
For full structure marks, the importing of any external modules (e.g. NumPy) should be done
at the top of your module. This should then be followed by all function definitions.
For full docstring marks, your docstrings should have the same style as discussed in lectures.
Important Marking Notes (read carefully):
• Please be aware that we will test your functions with an automated marker to check
that they work and produce the desired output(s), both with the data given in the question
and with different undisclosed data.
• .py files that include template functions, with correct inputs/outputs for all functions can
be found on Moodle.
• The template .py files also contain some simple tests of your functions. You are encouraged
to write your own tests of your functions in a separate file.
• To test that your code will pass through the automarker, download the files
test_student_code_Q*.py into the same folder as your solution .py files. Then run
test_student_code_Q*, which generates in your folder the file StudentCodeTestOutput_Q*.html.
Here, * is 1,2 or 3, depending on the question number. Open that file with a browser to see
exactly what tests were performed on your code. If your code fails to produce valid outputs
using this facility, then it will fail to run through the automarker, and you risk scoring low
output/plot marks.
• If your modules have filenames that differ to what has been explicitly asked for, then your
modules will not run through the automarker, and you risk scoring low output/plot marks.
Therefore, please do not add your username or student ID number to your filename.
• Every function you write should have a “docstring”, the contents of which could be
manually marked.
• Every function you write should have appropriate comments and a clear structure, which
could be manually marked.
• The automated tests for all questions will be ”timed out ” if they take too long to run. All
questions should only need at most a few seconds to run.
• When setting up a figure inside a function, you should first use the following syntax before
plotting anything:
fig = plt.figure()
fig can then be returned from the function as if it were any other variable.
MTHS2008
3 MTHS2008
Questions
1. For integration of (, ) over the rectangular region [, ]×[, ], we have already encountered
numerical quadrature of the form:
∫

∫


(, ) dd ≈
∫

[

∑
=0
̂(, )]
d
≈

∑
=0
[

∑
=0
̂(, )]
,
for some suitable evaluation points {}

=0 , {}

=0 and corresponding weights {}

=0 ,
{̂}

=0.
Suppose, instead, we would like to perform the integration
∫

∫
2()
1()
(, ) dd (1)
for some functions 1() and 2() that provide the limits in the -direction. We can apply a
numerical quadrature rule in both coordinate directions in this case, but notice that for each
quadrature point in the -direction, the quadrature points in the -direction will change.
Your task:
Complete the function found in multiple_integration.py (available on Moodle) with the
definition:
def double_integration(f,a,b,y1,y2,n)
• This function should find the approximation using the (composite) Simpson’s rule
with ≥ 2 strips in each coordinate direction.
• This function should take as input:
– A function f that itself will take two inputs x and y - the integrand;
– Floating points numbers a and b - the -limits of integration;
– Functions y1 and y2, that both take as input x - the -limits of integration;
– Integer n - the number of strips to use in each coordinate direction.
• The function as output:
– A floating point number integral - the (approximate) value of the integral (1).
Once you have written double_integration test it by running multiple_integration.py. You
should obtain the following:
integral1 = 3.5,
integral2 = 4.0.
Very tiny changes to these values are acceptable.
MTHS2008 Turn Over
4 MTHS2008
2. Suppose we wish to solve a system of initial value problems of the form:
du
d
() = u() + f(), 0 < ≤ , (2)
u(0) = u0, (3)
where u() is an unknown -dimensional vector function, is an × matrix and f() is a
known -dimensional vector function. Here u0 is the initial condition, which is a vector of
length .
The trapezoidal method for approximating solutions to (2)-(3) is given by
u+1 − u
Δ
=
1
2 (
(u + u+1) + f() + f(+1)) , = 0, 1, 2,… ,, (4)
where Δ = / is the time step, = Δ, for = 0, 1, 2,… , and u
≈ u().
Your task:
Complete the function found in time_steppers.py with the definition:
def trapezoidal_method(m,A,f,T,u0,N)
• This function should take as input:
– An integer m - the size of the ODE system;
– A sparse matrix A (type scipy.sparse.csr.csr_matrix) - the matrix from (2);
– A function f that itself takes an input and returns a vector of length - the function
f() from (2);
– A floating point number T - the final time;
– A vector u0 - the initial condition from (3);
– An integer N - the total number of time-steps to perform.
• The function should return as output:
– An array t of shape (N+1,), of all the time points , 0 ≤ ≤ .
– An array u of shape (N+1,m), where the nth row contains the approximation u found
using (4). Hence, the initial condition is stored in the 0th row.
• At each time step, a linear system solve should be performed using sparse.linalg.spsolve
from scipy.
• Hint: First rearrange (4) to make u+1 the subject and be very careful with what is
when doing u + u/2 = u, etc.
Once you have written trapezoidal_method.py test it by running time_steppers.py. You should
obtain the following:
t = [0. 0.2 0.4 0.6 0.8 1. ]
u =
[[0. 1. 0. 1. ]
[0.03216783 0.66853147 0.18601399 0.82237762]
[0.11636755 0.45521052 0.35305394 0.69189691]
[0.23320791 0.32862489 0.51450058 0.60991756]
[0.37140684 0.2676571 0.67671828 0.57296854]
[0.52420161 0.25749727 0.84244623 0.57574189]]
MTHS2008
5 MTHS2008
3. Suppose we wish to solve the following time dependent PDE for (, ):


−
2
2
+


+ = (), 0 ≤ ≤ , < < , (5)
where > 0, ≥ 0 and ≥ 0. The PDE is supplemented with initial condition
(, 0) = 0(), (6)
and boundary conditions
(, ) =


(, ) = 0. (7)
As usual, in the spatial dimension a mesh of + 1 equally spaced points should be used
so that
= + ℎ, = 0, 1,… ,,
with ℎ = ( − )/.
After finite difference semi-discretisation in space and letting () ≈ (, ), we obtain the
ODE system (2)-(3), with = + 1 and
=
⎛
⎜
⎜
⎜
⎜
⎝
1 0 0 0 … 0
0 … 0
0 … 0
⋮ ⋮ ⋱ ⋱ ⋱ …
0 … 0
0 … 0 0
⎞
⎟
⎟
⎟
⎟
⎠
, f() =
⎛
⎜
⎜
⎜
⎜
⎝
0
()
()
⋮
()
()
⎞
⎟
⎟
⎟
⎟
⎠
, (8)
where = /ℎ2 + /ℎ, = −(2/ℎ2 + /ℎ + ), = /ℎ2, = 2/ℎ2 + /ℎ and
u() =
⎛
⎜
⎜
⎜
⎜
⎜
⎝
0()
1()
2()
⋮
−1()
()
⎞
⎟
⎟
⎟
⎟
⎟
⎠
, u0 =
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
0(0)
0(1)
0(2)
⋮
0(−1)
0()
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
. (9)
In the time dimension, time-stepping scheme (4) is to be used with time-steps, so that
= Δ,
where Δ = /.
MTHS2008 Turn Over
6 MTHS2008
Your tasks:
(a) Complete the function found in finite_difference_1d.py with the definition:
def compute_A(a,b,c,Nx,h)
• This function should take as input:
– Floating point numbers a, b, c - the coefficients in (5);
– An integer Nx - the number of points in the -direction (minus 1).
– A floating point number h - the mesh width.
• The function should return:
– A sparse ( +1) × ( +1) matrix A (type scipy.sparse.csr.csr_matrix) - the
matrix in (8).
Once you have written compute_A.py, test it by running finite_difference_1d.py. You
should obtain the following:
A =
(0, 0) 1.0
(1, 0) 30.0
(1, 1) -56.0
(1, 2) 25.0
(2, 1) 30.0
(2, 2) -56.0
(2, 3) 25.0
(3, 2) 30.0
(3, 3) -56.0
(3, 4) 25.0
(4, 3) 30.0
(4, 4) -56.0
(4, 5) 25.0
(5, 4) 55.0
(5, 5) -56.0
As usual, small variations are permitted.
(b) Complete the function found in finite_difference_1d.py with the definition:
def pde_solver(a,b,c,xa,xb,T,Nx,Nt,f,u0)
• This function should take as input:
– Floating point numbers a, b, c - the coefficents in (5);
– Floating point numbers xa, xb - the domain limits;
– Floating point number T - the final time;
– An integer Nx - the number of points in the -direction (minus 1);
– An integer Nt - the number of time steps to perform;
– A function f, that itself takes input - the forcing function () from (5).
– A function u0, that itself takes input - the initial condition 0() from (6).
MTHS2008
7 MTHS2008
• The function should return:
– A numpy.ndarray u - an array of shape (Nt+1,Nx+1) that contains the approximate
solution to at all time points;
– A matplotlib.Figure.figure called fig - a surface plot of the approximate solution
to (, ).
• Your function should call compute_A from part (a) and trapezoidal_method from Q2.
• Hint: You will need to define a lambda function that implements f() from (8).
Once you have written pde_solver.py test it by running finite_difference_1d.py. You
should obtain the following:
u =
[[0. 0.30901699 0.58778525 0.80901699 0.95105652 1. ]
[0. 0.12848721 0.25445464 0.36592952 0.44834781 0.4838631 ]
[0. 0.11150322 0.19768906 0.26071885 0.2991047 0.31069468]
[0. 0.09130835 0.17014652 0.23104533 0.27154885 0.28863439]
[0. 0.09529189 0.16911275 0.22464435 0.26003658 0.27160544]]
Note: When marking your code, the model solution from Q2 will also be used, so even
if your Q2 is not working perfectly, you can obtain marks for Q3.
End of questions
MTHS2008 End