MATH2019 Introduction to Scientific Computation
— Coursework 2 (10%) —
Submission deadline: 3pm, Tuesday, 7 Dec 2021
Note: This is currently version 4 of the PDF document. (24th November, 2021)
Further questions will be added to this PDF in the upcoming weeks.
This coursework contributes 10% towards the overall grade for the module.
Rules:
• Each student is to submit their own coursework.
• You are allowed to work together and discuss in small groups (2 to 3 people), but you must write your
own coursework and program all code by yourself.
• Please be informed of the UoN Academic Misconduct Policy (incl. plagiarism, false authorship, and
collusion).
Coursework Aim and Coding Environment:
• In this coursework you will develop Python code related to algorithms that solve linear systems, and
you will study some of the algorithms’behaviour.
• You should write (and submit) plain Python code (.py), and 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).
How and Where to run Spyder:
• Spyder comes as part of the Anaconda package (recall that Jupyter is also part of Anaconda). Here
are three options on how to run Spyder:
(1) You can choose to install Anaconda on you personal device (if not done already).
(2) You can open Anaconda on any University of Nottingham computer.
(3) You can open a UoN Virtual Desktop on your own personal device, which is virtually the same as
logging onto a University of Nottingham computer, but through the virtual desktop. Then simply open
Anaconda. Here is further info on the UoN Virtual Desktop.
• The A18 Computer Room in the Mathematical Sciences building has a number of computers available,
as well as desks with dual monitors that you can plug into your own laptop.
• Getting started with Spyder? Open Spyder, and watch the basic intro video from the Spyder website
(click“Watch video”). You can also watch the recording (via Echo360) of Lecture 1, second hour, where
a brief demonstration was given.
Online 1-on-1 Support and Time-tabled Support Sessions (Thu 3-4pm, Wed 9-10am):
• You can work on the coursework whenever you prefer.
• We especially encourage you to work on it during the time-tabled computing and drop-in sessions
(Thursdays 3-4pm and Wednesdays 9-10am). You can visit the time-table allocated locations, or you
can obtain online 1-on-1 support via MS Teams, where PGR Student Teaching Assistants provide support.
• To obtain online 1-on-1 support during these sessions, please join the MATH2019 (21-22) MS Team
by using the code fnen7oo.
• For those of you that wish to go to the time-table allocated locations: Thursdays (3-4pm) is in Pharmacy
A06 (small seminar room with capacity of 33) and Wednesdays (9-10am) is in Clive Granger
(computer room with capacity of 68). You can find the lecturer there during these sessions.
School of Mathematical Sciences — Version: 24th November, 2021 — Lecturer: Dr Kris van der Zee Page 1 of 11
Piazza:
• You are allowed and certainly encouraged(!) to also ask questions using Piazza to obtain clarification
of the coursework questions, as well as general Python queries.
• However, when using Piazza, please ensure that you are not revealing any answers to others. Also,
please ensure that you are not directly revealing any code that you wrote. Doing so is considered
Academic Misconduct.
• When in doubt, please simply attend a Support Sessions (Thursdays 3-4pm, Wednesdays 9-10am).
Some Further Advice:
• You are expected to have basic familiarity with Python, in particular: logic and loops, functions, NumPy
and matplotlib.pyplot. Note that it will always be assumed that the package numpy is imported as np,
and matplotlib.pyplot as plt.
• Helpful resources (web): Python 3 Online Documentation – Spyder IDE (integrated development
editor) – NumPy User Guide – Matplotlib Usage Guide
• Helpful resources (Core Programming Course): Moodle Page Core Programming 2020-2021 –
Moodle Page Core Programming 2021-2022 – Content of Core Programming (20-21): See MATH2019
Moodle Page for Jupyter Notebooks A–J [MATH2019 uses: A, B, C(§1.4), E(§2,§4), F, G]
• Write your code as neatly and read-able as possible so that it is easy to follow. Add some comments
to your code that indicate what a piece of code does. Frequently run your code to check for (and
immediately resolve) mistakes and bugs.
• Coursework Questions with a“?”are more tedious, but can be safely skipped, as they don’t affect
follow-up questions.
Submission Procedure:
• Submission will open after 25 November 2021.
• To submit, simply upload the requested .py-files on Moodle. (Your submission will be checked for
plagiarism using turnitin.)
• Your work will be marked (mostly) automatically: This is done by running your functions and comparing
their output against the true results.
Getting Started:
• Download the contents of the“Coursework 2 Pack”folder from Moodle into a single folder.
(More files may be added in later weeks.)
I Gaussian elimination without pivoting
Let A ∈ R
n×n and b ∈ R
n
. Suppose that A is such that there exists a unique solution
x ∈ R
n
to Ax = b and that forward elimination can be performed without any
row interchanges to arrive at an augmented matrix to which backward substitution
can be applied to find x.
1 Open the file systemsolvers.py. Note that this file already contains an unfinished
function with the following definition:
def no_pivoting (A ,b ,n , c)
This function returns the output M which is of type numpy.ndarray and has shape
(n, n + 1). The input A is of type numpy.ndarray and has shape (n, n) and represents
the square matrix A. The input b is of type numpy.ndarray and has
Coursework 2 Page 2 of 11
shape (n, 1) and represents the column vector b. The input n is a positive integer
an integer ≥ 2,
← Sentence
corrected on
4 Nov 10am
The input c is an integer such that 1 ≤ c ≤ n − 1, and it is used
to (prematurely) stop the forward elimination algorithm; see details below.
• A pseudocode algorithm for performing forward elimination without row interchanges
on the matrix M ∈ R
n×n+1 is the following:
For i = 1 to n − 1 do
For j = i + 1 to n do
Set g = mj,i/mi,i
Set mj,i = 0
For k = i + 1 to n + 1 do
Set mj,k = mj,k − gmi,k
End do
End do
End do
In the above pseudocode algorithm, the entries of the matrix M are mi,j for
i = 1, 2, . . . , n and j = 1, 2, . . . , n + 1. However, the elements in M (a
numpy.ndarray with shape (n, n + 1)) are M[i,j] for i = 0, 1, . . . , n − 1 and
j = 0, 1, . . . , n.
• Complete the no pivoting function so that the output M is a 2-D array representing
the augmented matrix arrived at by starting from the augmented matrix
[A b] and performing forward elimination without row interchanges until all of
the entries below the main diagonal in the first c columns are 0.
• The backward.py file contains a finished function with the following definition:
def backward_substitution (M , n)
The input M is of type numpy.ndarray and has shape (n, n+1) and represents an
augmented matrix of the form [U v] where U ∈ R
n×n
is an upper triangular
matrix with all entries on its main diagonal being nonzero and v ∈ R
n
. The input
n is a positive integer. The output x is of type numpy.ndarray and has shape
(n, 1), and represents the solution x to Ux = v.
• Add in your file systemsolvers.py the following function definition:
def no_pivoting_solve (A ,b , n)
The input A is of type numpy.ndarray and has shape (n, n) and represents the
matrix A. The input b is of type numpy.ndarray and has shape (n, 1) and represents
the vector b. The input n is a positive integer.
• Complete the no pivoting solve function so that the output x is of type
numpy.ndarray and has shape (n, 1), and represents the solution x to Ax = b
computed using Gaussian elimination without pivoting. ← Sentence
extended on
11 Nov 11am
The no pivoting solve
Coursework 2 Page 3 of 11
function can call the no pivoting and backward substitution functions in order
to do this.
• Also add a brief description at the top of your definition, the so-called
doc string. This description should become visible whenever one types:
help(ss.no pivoting solve) or ss.no pivoting solve?.
• Test your no pivoting and no pivoting solve functions by running the
main.py file, which contains the following:
%% Question 1
Initialise
A = np . array ([[1 , -5 ,1] ,[10 ,0.0 ,20] ,[5 ,10 , -1]])
b = np . array ([[7] ,[6] ,[4]])
n = 3;
Run no_pivoting
M = ss . no_pivoting (A ,b ,n ,1)
Print output
print (M)
Run no_pivoting
M = ss . no_pivoting (A ,b ,n ,2)
Print output
print (M)
Solve Ax=b
x = ss . no_pivoting_solve (A ,b , n)
Print output
print (x)
Help
help (ss . no_pivoting_solve)
For the above A, b and n, with c=1, the output from the no pivoting function
should be:
[[1. -5. 1. 7.]
[0. 50. 10. -64.]
[0. 35. -6. -31.]]
For the above A, b and n, with c=2, the output from the no pivoting function
should be:
[[1. -5. 1. 7.]
[0. 50. 10. -64.]
[0. 0. -13. 13.8]]
(The below has been added on 11 Nov 2021.)
I
Marks can be obtained for your no pivoting function definition for correctly gener- [13 / 40]
ating the required output, for certain set(s) of inputs for {A, b, n, c}. The correctness
of the following will be checked:
• The type of output M
Coursework 2 Page 4 of 11
• The np.shape of output M
• The values of output M
Marks can be obtained for your no pivoting solve function definition for correctly
generating the required output, for certain set(s) of inputs for {A, b, n}. The
correctness of the following will be checked:
• The type of output x
• The np.shape of output x
• The values of output x
• The output of“help(ss.no pivoting solve)”.
To test your code, run test student code which generates in your folder the file
StudentCodeTestOutput.html. Open that file with a web browser to see exactly
what tests were performed on your code.
Note that in marking your work, different input(s) may be used.
I Gaussian elimination with partial pivoting
Let A ∈ R
n×n and b ∈ R
n
. Suppose that A is such that there exists a unique
solution x ∈ R
n
to Ax = b and that forward elimination with partial pivoting can be
performed to arrive at an augmented matrix to which backward substitution can be
applied to find x.
2 A pseudocode algorithm for finding the smallest positive integer p such that
|vp| = max
j∈{1,2,…,n}
|vj
|
is the following:
Set m = 0
Set p = 1
For j = 1 to n do
If |vj
| > m
Set m = |vj
|
Set p = j
End If
End do
• Add in your file systemsolvers.py the following function definition:
def find_max (M ,n , i)
The input M is of type numpy.ndarray and has shape (n, n + 1). The input n
is an integer such that n ≥ 2. The input i is a nonnegative integer such that
i ≤ n − 2.
• Complete the find max function so that the output p is of type int and is the
Coursework 2 Page 5 of 11
smallest integer such that p ≥ i and
|M[p, i]| = max
j∈{i,i+1,…,n−1}
|M[j, i]|.
• Add in your file systemsolvers.py the following function definition:
def partial_pivoting (A ,b ,n , c)
The input A is of type numpy.ndarray and has shape (n, n) and represents the
square matrix A. The input b is of type numpy.ndarray and has shape (n, 1)
and represents the column vector b. The input n is an integer such that n ≥ 2.
The input c is an integer such that 1 ≤ c ≤ n − 1, and it is used to (prematurely)
stop the forward elimination with partial pivoting algorithm; see details below.
• Complete the partial pivoting function so that the output M is of type
numpy.ndarray and has shape (n, n + 1), and represents the augmented matrix
arrived at by starting from the augmented matrix [A b] and performing
forward elimination with partial pivoting (as described in Lecture 5) until
all of the entries below the main diagonal in the first c columns are 0. The
partial pivoting function can call the find max function in order to do this.
Note that M[[i,p],:]=M[[p,i],:] can be used to interchange the row M[i,:]
and the row M[p,:].
• Add in your file systemsolvers.py the following function definition:
def partial_pivoting_solve (A ,b , n)
The input A is of type numpy.ndarray and has shape (n, n) and represents the
matrix A. The input b is of type numpy.ndarray and has shape (n, 1) and represents
the vector b. The input n is an integer such that n ≥ 2.
• Complete the partial pivoting solve function so that the output x is
of type numpy.ndarray and has shape (n, 1), and represents the solution
x to Ax = b computed using Gaussian elimination with partial pivoting.
The partial pivoting solve function can call the partial pivoting and
backward substitution functions in order to do this.
• Test your find max, partial pivoting and partial pivoting solve functions
by running the main.py file, which has been updated to contain the following:
%% Question 2
import systemsolvers as ss
Initialise
A = np . array ([[1 , -5 ,1] ,[10 ,0.0 ,20] ,[5 ,10 , -1]])
b = np . array ([[7] ,[6] ,[4]])
n = 3
Coursework 2 Page 6 of 11
Run find_max
p = ss . find_max (np . hstack (( A , b) ) ,n ,0)
Print output
print (p)
Run partial_pivoting
M = ss . partial_pivoting (A ,b ,n ,1)
Print output
print (M)
Run find_max
p = ss . find_max (M ,n ,1)
Print output
print (p)
Run partial_pivoting
M = ss . partial_pivoting (A ,b ,n ,2)
Print output
print (M)
Solve Ax=b
x = ss . partial_pivoting_solve (A ,b , n)
Print output
print (x)
For the above A, b and n, with c=1, the output from the partial pivoting
function should be:
[[10. 0. 20. 6.]
[0. -5. -1. 6.4]
[0. 10. -11. 1.]]
I
Marks can be obtained for your find max function definition for correctly generat- [13 / 40]
ing the required output, for certain set(s) of inputs for {M, n, i}. The correctness of
the following will be checked:
• The type of output p
• The value of output p
Marks can be obtained for your partial pivoting function definition for correctly
generating the required output, for certain set(s) of inputs for {A, b, n, c}. The correctness
of the following will be checked:
• The type of output M
• The np.shape of output M
• The values of output M
Marks can be obtained for your partial pivoting solve function definition for
correctly generating the required output, for certain set(s) of inputs for {A, b, n}. The
correctness of the following will be checked:
Coursework 2 Page 7 of 11
• The values of output x.
To test your code, run test student code which has been updated to include tests
for Q2 and generates in your folder the file StudentCodeTestOutput.html. Open
that file with a web browser to see exactly what tests were performed on your code.
Note that in marking your work, different input(s) may be used.
(The below has been added on 18 Nov 2021.)
I Doolittle factorisation (Note: Question 3 is not required in order to do Question 4.)
Let A ∈ R
n×n be such that Theorem 6.19 in [Burden, Faires & Burden, Numerical
Analysis, 10E] (see also Lecture 6 Slides) can be applied to obtain a lower triangular
matrix L ∈ R
n×n with the entries on its main diagonal all being 1 and an upper
triangular matrix U ∈ R
n×n
such that A = LU.
3? Add in your file systemsolvers.py the following function definition:
def Doolittle (A , n)
The input A is of type numpy.ndarray and has shape (n, n) and represents the
matrix A. The input n is an integer such that n ≥ 2.
• Complete the Doolittle function so that
return L, U
is used to output both L and U. The output L is of type numpy.ndarray and
has shape (n, n), and represents the lower triangular matrix L in the Doolittle
factorisation of A. The output U is of type numpy.ndarray and has shape (n, n),
and represents the upper triangular matrix U in the Doolittle factorisation of A.
See np.identity? for a useful starting point for L.
Note: You are not allowed to use an existing Python function that obtains a
factorization of A. Hence, there is no need to import a package such as scipy.
Instead, you are supposed to directly implement the result for L and U stated in
Theorem 6.19.
• Test your Doolittle function by running the main.py file, which has been
updated to contain the following:
%% Question 3
import systemsolvers as ss
Initialise
A = np . array ([[1 ,1 ,0.0] ,[2 ,1 , -1] ,[0 , -1 , -1]])
n = 3
Run Doolittle
L , U = ss . Doolittle (A , n)
Print output
print (L)
print (U)
Coursework 2 Page 8 of 11
I
Marks can be obtained for your Doolittle function definition for correctly gener- [6 / 40]
ating the required output, for certain set(s) of inputs for {A, n}. The correctness of
the following will be checked:
• The type of output L and output U
• The np.shape of output L and output U
• The values of output L and output U.
To test your code, run test student code which has been updated to include tests
for Q3 and generates in your folder the file StudentCodeTestOutput.html. Open
that file with a web browser to see exactly what tests were performed on your code.
Note that in marking your work, different input(s) may be used.
(The below has been added on 24 Nov 2021.)
I Gauss–Seidel method
Let A ∈ R
n×n and b ∈ R
n
. Suppose that A is such that there exists a unique
solution x ∈ R
n
to Ax = b and that all of the entries on the main diagonal of A are
nonzero. Let x
(k) be the approximation to x obtained after performing k iterations of
the Gauss–Seidel method, starting from the initial approximation x
(0)
.
4? Add in your file systemsolvers.py the following function definition:
def Gauss_Seidel (A ,b ,n , x0 , tol , maxits)
The input A is of type numpy.ndarray and has shape (n, n) and represents the
matrix A. The input b is of type numpy.ndarray and has shape (n, 1) and represents
the vector b. The input n is an integer such that n ≥ 2. The input x0 is of
type numpy.ndarray and has shape (n, 1) and represents the initial approximation
x
(0). The input tol is a positive real number. The input maxits is a positive
integer.
• Complete the Gauss Seidel function so that the output x is of type
numpy.ndarray and has shape (n, 1), and represents the approximation x
(k)
where k is the smallest positive integer such that k ≤ maxits and
kb − Ax(k)
k∞ < tol
if such a k exists. Note that np.max(np.abs(b-np.matmul(A,x))) can be used
to compute kb − Ax(k)k∞. If no such k exists then the output x should be of
type str and should be the string “Desired tolerance not reached after
maxits iterations have been performed.”.
Hint: You can choose to implement directly the element-based form of the
Gauss-Seidel method (see Lecture 7 Slides and Notes), or you can choose to
implement a vector-based form of the Gauss-Seidel method (see also Lecture
7), in which case the following functions could be useful: np.tril, np.triu,
Coursework 2 Page 9 of 11
np.diag, np.linalg.inv, np.matmul or the’@’operator for matrix multiplication.
• Test your Gauss Seidel function by running the main.py file, which has been
updated to contain the following:
%% Question 4
import systemsolvers as ss
Initialise
A = np . array ([[4 , -1 ,0] ,[-1 ,8 , -1] ,[0 , -1 ,4]])
b = np . array ([[48] ,[12] ,[24]])
n = 3
x0 = np . array ([[1.0] ,[1] ,[1]])
tol = 1e -2
Run Gauss_Seidel
x = ss . Gauss_Seidel (A ,b ,n , x0 , tol ,3)
Print output
print (x)
Run Gauss_Seidel
x = ss . Gauss_Seidel (A ,b ,n , x0 , tol ,4)
Print output
print (x)
For the above A, b, n, x0 and tol, with maxits=3, the output from the
Gauss Seidel function should be:
Desired tolerance not reached after maxits iterations have
been performed .
For the above A, b, n, x0 and tol, with maxits=4, the output from the
Gauss Seidel function should be:
[[12.99917603]
[3.99979401]
[6.9999485]]
I
Marks can be obtained for your Gauss Seidel function definition for correctly gen- [8 / 40]
erating the required output, for certain set(s) of inputs for {A, b, n, x0, tol, maxits}.
The correctness of the following will be checked:
• The type of output x
• The np.shape of output x when applicable
• The values of output x
To test your code, run test student code which has been updated to include tests
for Q4 and generates in your folder the file StudentCodeTestOutput.html. Open
Coursework 2 Page 10 of 11
that file with a web browser to see exactly what tests were performed on your code.
Note that in marking your work, different input(s) may be used.
I Finished!
Submission Procedure:
• Submission will open after 25 November 2021.
• To submit, simply upload your file systemsolvers.py on Moodle using the“CW2
Submission”Assignment. (Your submission will be checked for plagiarism using
turnitin.)
• Do not change the name of the systemsolvers.py file. The name of the file that you
upload must be systemsolvers.py
Please don’t forget to click the submit button. Thank you!
Coursework 2 Page 11 of 11