代写COMP2421编程、代做data编程、Python程序语言调试

日期：2020-11-10 11:22

School of Computing, University of Leeds, 2020–2021

COMP2421: Numerical Computation

Coursework Two: 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 (typed or hand-written if fine) should be submitted electronically via

Minerva by 10:00 on Monday 8th November 2020.

Additional instructions: Please include code snipets (rather than entire files) as appropriate

as part of your answers. Please upload your submission as a single PDF file.

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 = n

2

) 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).