SWJTU-Leeds Joint School, 2019–2020XJCO2421: Numerical ComputationCoursework Five: Large Linear Systems – Python Version (Summative)These exercises are intended to increase your understanding of the application and limitationsof GE and iterative methods in computational modelling. Both questions are assessedand your solutions (hand-written if fine) should be submitted electronically via Minerva by16:00 on Monday 11th November 2019.Q1. For this question you are required to modify the functions gauss elimination andupper triangular solve in the file matrixSolve.py that were introduced in Lectures 5and 6, in order to make them return the number of arithmetic operations that each of themhave executed (in addition to the data that they currently return). You should then writea new function called gauss elim count which has a single integer parameter, n say, asits argument and returns two integers: the number of operations required for the forwardelimination of an n × n system, and the number of operations required for the backwardsubstitution. Finally, you should provide an analysis of how the work grows as the value ofn 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 thatthe diagonal is nonzero, i.e. A(i,i) \= 0.ARGUMENTS: A n x n matrixb right hand side column n-vectorRETURNS: A upper triangular n x n matrixb modified column n-vectorcount count of the number of operations in the inner loopNote that you will need to insert your counter of the number of operations in theinner-most loop of the elimination. Furthermore, for this exercise you may regard amultiplication, subtraction and assignment as a single operation (e.g. x = x − y ∗ zcan 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.e. A(i,i) \= 0.ARGUMENTS: A upper triangular n x n matrixb right hand side column n-vectorRETURNS: x column n-vector solutioncount number of operations in the inner loopNote that, once again, you will need to insert a counter for the number of operationsin 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 Eliminationand solving the resulting upper triangular system, but return the numberof operations executed by the Elimination and the backward substitution.ARGUMENTS: n dimension of systemRETURNS: count1 operations in forward elimination (GE)count2 operations in backward substitutionThis should make use of numpy’s rand function to generate a random n × n matrixand a random right-hand side vector, and then use your modified versions ofgauss elimination and upper triangular solve to solve the system and evaluatecount1 and count2. [4 marks](d) Use your function gauss elim count to produce a table that displays the forwardelimination 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 andshow the ratio of the operation count for the current value of n divided by the countfor the previous value of n. [4 marks](e) In each case (forward elimination and backward substitution) what do your results tellyou about how the number of operations grows each time n is doubled? What do youdeduce 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 forparts (a) to (c), your computed results for part (d) and your discussion of these results forpart (e).Q2. The Python script file heat.py can be used to model and approximate the temperaturedistribution in a two-dimensional square flat plate. This is achieved by approximating theplate as a grid of m × m connected points as illustrated in Fig. 1. The temperature at eachpoint in the grid is to be stored in the array T emp: with the values on the boundary beingspecified and the values in the interior being computed based upon a model which statesthat 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 isstored in the variable n2 in the program. The n2×n2 array A is used to store the coefficientmatrix for the resulting system of linear equations, and the vector b is used to store theright-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 unknoXJCO2421代做、Python程序设计调试、代写Linewn temperatures the programsolves 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 arethen copied into the corresponding entries of the array T emp so that this array now storesthe approximate temperatures throughout the whole of the plate. The output is then asurface plot of the overall temperature distribution along with the approximated value atthe 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 sidevector b are influenced by the known temperature on the boundary. Finally, explainwhy no row of the array A can ever have more than 5 non-zero entries (regardless ofthe choice of m). [5 marks](b) Run the script heat.py with the given temperature distribution on the boundary forgradually 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). Producea table of results that contains the following data: the choice of m, the number ofunknowns in the resulting linear system, the estimated temperature at the centre ofthe square, the error in the temperature at the centre of the square (for the givenboundary values you may assume that the exact answer should be 19.8671) and theexecution 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 ofm and how this, in turn, affects the execution time taken (NB the solution time form = 17 is so small that you should ignore this and just consider the other three caseswhen looking for a pattern!). [7 marks](c) Given that Python stores real numbers (float type) in 64 bits (i.e. 8 bytes) approximatelyhow much memory is required for the array A in the cases m = 129, and howmuch memory would be required to store A when m = 257? Conversely, given thatthere are at most 5 non-zero entries in each row of A (see part (a) above) and thatonly these non-zero entries need be stored (in 128 bits each: 64 bits for the entry and32 bits each for the row and column number), how much memory would be requiredto 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 aniterative solver based upon Jacobi or Gauss-Seidel iteration (using the form that waspresented in Lecture 9: where a convergence tolerance is used) - remember that youwill 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 tothe given tolerance (10−5) with each of the iterative methods (Jacobi and Gauss-Seideliteration) for systems of equations with m = 9, m = 17 and m = 33. What canyou say about the relative number of iterations in each case and about the way inwhich the iteration count grows with m? If you decrease (or increase) the convergencetolerance 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. Thisis obtained by settingu(k+1)j = wuj + (1 − w)u(k)jwhere ujis the value for u(k+1)jthat would have been obtained from Gauss-Seideliteration and 0 is the same as Gauss-Seidel iteration however when w > 1 it generally convergesfaster 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 newfunction called sor new that has the following functionality:Solve the system A u = b using SOR iterationARGUMENTS: A k x k matrixu k-vector storing initial estimateb k-vector storing right-hand sidek integer dimension of systemtol real number providing the required convergence tolerancew weight factor for the iterationRESULTS: u k-vector storing solution[5 marks](g) Now use your new function sor new to solve the heat problem for m = 33 with aconvergence tolerance of 10−5, trying different values of w (between 1.0 and 2.0). Statehow many iterations are required to reach a converged solution for each value of w thatyou try and see if you can find the best possible choice of value for this problem? Nowconsider the problem using a different grid (e.g. m = 17): what is the optimal choiceof w in this case? Has it changed? [5 marks]In addition to your written answers to each part of this question please hand in commentedlistings of your modified code for parts (b) and (d) (you can hand in a single listing with themodifications for both parts: the timing and the call to one of the iterative solvers), as wellas a a commented listing of your SOR function for part (f).转自：http://www.3daixie.com/contents/11/3444.html