SWJTU-Leeds Joint School, 2019–2020
XJCO2421: Numerical Computation
Coursework Five: Large Linear Systems – Python Version (Summative)
These exercises are intended to increase your understanding of the application and limitations
of GE and iterative methods in computational modelling. Both questions are assessed
and your solutions (hand-written if fine) should be submitted electronically via Minerva by
16:00 on Monday 11th November 2019.
Q1. For this question you are required to modify the functions gauss elimination and
upper triangular solve in the file matrixSolve.py that were introduced in Lectures 5
and 6, in order to make them return the number of arithmetic operations that each of them
have executed (in addition to the data that they currently return). You should then write
a new function called gauss elim count which has a single integer parameter, n say, as
its argument and returns two integers: the number of operations required for the forward
elimination of an n × n system, and the number of operations required for the backward
substitution. Finally, you should provide an analysis of how the work grows as the value of
n is increased. More specific instructions are as follows.
(a) Modify gauss elimination so that it has the following functionality:
"""
Reduce the system A x = b to upper triangular form, assuming that
the diagonal is nonzero, i.e. A(i,i) \= 0.
ARGUMENTS: A n x n matrix
b right hand side column n-vector
RETURNS: A upper triangular n x n matrix
b modified column n-vector
count count of the number of operations in the inner loop
"""
Note that you will need to insert your counter of the number of operations in the
inner-most loop of the elimination. Furthermore, for this exercise you may regard a
multiplication, subtraction and assignment as a single operation (e.g. x = x − y ∗ z
can be treated as a single operation). [4 marks]
(b) Modify upper triangular solve so that it has the following functionality:
"""
Solve the system A x = b where A is assumed to be upper triangular,
i.e. A(i,j) = 0 for j < i, and the diagonal is assumed to be nonzero,
i.e. A(i,i) \= 0.
ARGUMENTS: A upper triangular n x n matrix
b right hand side column n-vector
RETURNS: x column n-vector solution
count number of operations in the inner loop
"""
Note that, once again, you will need to insert a counter for the number of operations
in the inner-most loop. [4 marks]
(c) Write a new function called gauss elim count with the following functionality:
"""
Solve a nxn example system A x = b by first using Gaussian Elimination
and solving the resulting upper triangular system, but return the number
of operations executed by the Elimination and the backward substitution.
ARGUMENTS: n dimension of system
RETURNS: count1 operations in forward elimination (GE)
count2 operations in backward substitution
"""
This should make use of numpy’s rand function to generate a random n × n matrix
and a random right-hand side vector, and then use your modified versions of
gauss elimination and upper triangular solve to solve the system and evaluate
count1 and count2. [4 marks]
(d) Use your function gauss elim count to produce a table that displays the forward
elimination operation counts and the backward substitution operation counts for n =
[2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]. In each case (except n = 2) also calculate and
show the ratio of the operation count for the current value of n divided by the count
for the previous value of n. [4 marks]
(e) In each case (forward elimination and backward substitution) what do your results tell
you about how the number of operations grows each time n is doubled? What do you
deduce from this about the computational complexity of each of these algorithms?
[4 marks]
Please hand in commented listings of the three functions that you have modified/written for
parts (a) to (c), your computed results for part (d) and your discussion of these results for
part (e).
Q2. The Python script file heat.py can be used to model and approximate the temperature
distribution in a two-dimensional square flat plate. This is achieved by approximating the
plate as a grid of m × m connected points as illustrated in Fig. 1. The temperature at each
point in the grid is to be stored in the array T emp: with the values on the boundary being
specified and the values in the interior being computed based upon a model which states
that the temperature at each point is the average of the temperature at its four neighbours.
Figure 1: Illustration of the grid model for the temperature distribution in a square plate:
in this example m = 9, the number of unknown temperatures in each direction (n = m−2)
is 7 and the total number of unknowns (n2 = n2) is 49.
The total number of unknown temperature values is therefore (m −2)×(m −2) which is
stored in the variable n2 in the program. The n2×n2 array A is used to store the coefficient
matrix for the resulting system of linear equations, and the vector b is used to store the
right-hand side of the system. The unknowns are numbered from top left to bottom right,
as illustrated in Fig. 1.
Having created the system of n2 equations for the n2 unknown temperatures the program
solves this system using Python’s built-in implementation of LU factorisation: u =
np.linalg.solve(A,b). The result is initially stored in the vector u but these values are
then copied into the corresponding entries of the array T emp so that this array now stores
the approximate temperatures throughout the whole of the plate. The output is then a
surface plot of the overall temperature distribution along with the approximated value at
the middle of the plate (assuming that m is odd so that there is a grid point at the centre).
(a) Explain why the first row of the array A always has 3 entries which are non-zero,
whereas the second row of the array A always has 4 non-zero entries (provided m ≥ 5).
In each of these cases explain how the first and second entries of the right-hand side
vector b are influenced by the known temperature on the boundary. Finally, explain
why no row of the array A can ever have more than 5 non-zero entries (regardless of
the choice of m). [5 marks]
(b) Run the script heat.py with the given temperature distribution on the boundary for
gradually increasing values of m (specifically, for m = 17, m = 33, m = 65 and m = 129
– be patient when executing this last case! – always using odd values for m). Produce
a table of results that contains the following data: the choice of m, the number of
unknowns in the resulting linear system, the estimated temperature at the centre of
the square, the error in the temperature at the centre of the square (for the given
boundary values you may assume that the exact answer should be 19.8671) and the
execution time to solve the problem (you should make use of time.perf counter()
in Python 3). Discuss how the accuracy of the solution depends upon the choice of
m and how this, in turn, affects the execution time taken (NB the solution time for
m = 17 is so small that you should ignore this and just consider the other three cases
when looking for a pattern!). [7 marks]
(c) Given that Python stores real numbers (float type) in 64 bits (i.e. 8 bytes) approximately
how much memory is required for the array A in the cases m = 129, and how
much memory would be required to store A when m = 257? Conversely, given that
there are at most 5 non-zero entries in each row of A (see part (a) above) and that
only these non-zero entries need be stored (in 128 bits each: 64 bits for the entry and
32 bits each for the row and column number), how much memory would be required
to store A in sparse format when m = 129 or m = 257 respectively? [8 marks]
(d) Now replace the built-in Python solver (“u= np.linalg.solve(A,b)”) by a call to an
iterative solver based upon Jacobi or Gauss-Seidel iteration (using the form that was
presented in Lecture 9: where a convergence tolerance is used) - remember that you
will now require an initial guess to the solution (e.g. the zero vector). [3 marks]
(e) Produce a table of results to show how many iterations are required to converge to
the given tolerance (10−5
) with each of the iterative methods (Jacobi and Gauss-Seidel
iteration) for systems of equations with m = 9, m = 17 and m = 33. What can
you say about the relative number of iterations in each case and about the way in
which the iteration count grows with m? If you decrease (or increase) the convergence
tolerance by successive factors of ten what happens to the iteration count in each case?
[7 marks]
(f) There is a modification of the Gauss-Seidel scheme that can converge faster still. This
is obtained by setting
u(k+1)j = wuj + (1 − w)u(k)j
where uj
is the value for u(k+1)j
that would have been obtained from Gauss-Seidel
iteration and 0 < w < 2 is known as a relaxation parameter. When w = 1 the method
is the same as Gauss-Seidel iteration however when w > 1 it generally converges
faster than Gauss-Seidel. For this reason the iteration is known as Successive OverRelaxation,
or SOR for short. Modify the function gauss seidel new to create a new
function called sor new that has the following functionality:
"""
Solve the system A u = b using SOR iteration
ARGUMENTS: A k x k matrix
u k-vector storing initial estimate
b k-vector storing right-hand side
k integer dimension of system
tol real number providing the required convergence tolerance
w weight factor for the iteration
RESULTS: u k-vector storing solution
"""
[5 marks]
(g) Now use your new function sor new to solve the heat problem for m = 33 with a
convergence tolerance of 10−5
, trying different values of w (between 1.0 and 2.0). State
how many iterations are required to reach a converged solution for each value of w that
you try and see if you can find the best possible choice of value for this problem? Now
consider the problem using a different grid (e.g. m = 17): what is the optimal choice
of w in this case? Has it changed? [5 marks]
In addition to your written answers to each part of this question please hand in commented
listings of your modified code for parts (b) and (d) (you can hand in a single listing with the
modifications for both parts: the timing and the call to one of the iterative solvers), as well
as a a commented listing of your SOR function for part (f).