MATH4063代写、C++编程设计代做

日期：2023-11-08 07:43

1 MATH4063

The University of Nottingham

SCHOOL OF MATHEMATICAL SCIENCES

AUTUMN SEMESTER 2022-2023

MATH4063 - SCIENTIFIC COMPUTING AND C++

Coursework 1 - Released 30th October 2023, 4pm

Your work should be submitted electronically via the MATH4063 Moodle page by 12noon on Monday 20th

November (unless you have arranged an extension). Since this work is assessed, your submission must be

entirely your own work (see the University’s policy on Academic Misconduct). Submissions up to five working

days late will be marked, but subject to a penalty of 5% of the maximum mark per working day.

The marks for each question are given by means of a figure enclosed by square brackets, eg [20]. There are

a total of 100 marks available for the coursework and it contributes 45% to the module. The marking rubric

available on Moodle will be applied to each full question to further break down this mark.

You are free to name the functions you write as you wish, but bear in mind these names should be meaningful.

Functions should be grouped together in .cpp files and accessed in other files using correspondingly named

.hpp files.

All calculations should be done in double precision.

A single zip file containing your full solution should be submitted on Moodle. This zip file should contain three

folders called main, source and include, with the following files in them:

main:

? q1d.cpp

? q2c.cpp

? q3b.cpp

? q4b.cpp

source:

? vector.cpp

? dense_matrix.cpp

? csr_matrix.cpp

? linear_algebra.cpp

? finite_volume.cpp

include:

? vector.hpp

? dense_matrix.hpp

? csr_matrix.hpp

? linear_algebra.hpp

? finite_difference.hpp

Prior to starting the coursework, please download the CW1_templates.zip from Moodle and extract the files.

More information about the contents of the files included in the template is given in the questions below.

Hint: When using a C++ struct with header files, the whole struct needs to be defined fully in the header file,

and the header file included in the corresponding .cpp file. Include guards should also be used.

MATH4063 Turn Over

2 MATH4063

In this coursework you will build a 2D finite volume solver for the following PDE boundary value problem

????? + ? ? (b??) = ?? (??, ??) ∈ ??, (1)

?? = ??, (??, ??) ∈ ????, (2)

where ?? ∶ ?? → ?, ?? ∶ ???? → ? and b ∶ ?? → ?2

.

In order to solve this problem, you will first define a sparse matrix structure, then write functions to apply

the GMRES linear algebra solver and finally build and solve the linear system arising from the finite volume

approximation of (1)-(2).

1. Matrices arising from the discretisation of partial differential equations using, for example, finite volume

methods, are generally sparse in the sense that they have many more zero entries than nonzero ones.

We would like to avoid storing the zero entries and only store the nonzero ones.

A commonly employed sparse matrix storage format is the Compressed Sparse Row (CSR) format. Here,

the nonzero entries of an ?? × ?? matrix are stored in a vector matrix_entries, the vector column_no gives

the column position of the corresponding entries in matrix_entries, while the vector row_start of length

??+1 is the list of indices which indicates where each row starts in matrix_entries. For example, consider

the following:

?? =

?

?

?

?

8 0 0 2

0 3 1 0

0 0 4 0

6 0 0 7

?

?

?

?

?

matrix_entries = (8 2 3 1 4 6 7)

column_no = (0 3 1 2 2 0 3)

row_start = (0 2 4 5 7)

Note, in the above, C++ indexing has been assumed, i.e, indices begin at 0.

(a) In csr_matrix.hpp, define a C++ struct called csr_matrix to store a matrix in CSR format. In

addition to matrix_entries, column_no and row_start, you should store the number of rows of the

matrix explicitly.

(b) In csr_matrix.cpp, write a C++ function that will set up the matrix ?? from above in CSR format.

Remember, if you are using dynamically allocated memory, then you should also have corresponding

functions that will deallocate the memory you have set up.

(c) In csr_matrix.cpp, write a C++ function that takes as input a matrix ?? stored in CSR format and a

vector x and computes the product ??x. The prototype for your function should be:

void MultiplyMatrixVector ( csr_matrix & matrix ,double* vector ,

double* productVector )

Hence, the input vector and the output productVector should be pointers to dynamically allocated

arrays. In particular, it should be assumed that productVector has been preallocated to the correct

size already.

(d) By setting a vector x = (4, ?1, 3, 6)?, write a test program in q1d.cpp to compute and print to the

screen the product ??x, where ?? is the matrix given above.

[20 marks]

MATH4063

3 MATH4063

2. Suppose we wish to find x ∈ ???

such that

??x = b, (3)

where ?? is an ?? × ?? matrix and b ∈ ???

.

One algorithm for solving this problem is the (restarted) Generalised Minimal RESidual (GMRES) algorithm.

The method is too complicated to explain here, but works to quickly find approximations x?? = x0 + y??

where y?? ∈ ???? ∶= Span{??q0

, ??2q0 … ????q0

} for ?? = 1, 2, …. y?? is chosen to minimise the residual

‖b ? ??x??‖2

.

Here x0

is some initial guess vector and q0

is the normed initial residual

q0 =

b ? ??x0

‖b ? ??x0‖2

.

???? is called a Krylov subspace of ??.

The algorithm stops when ‖b ? ??x??‖2 < tol for some termination tolerance tol. As the method becomes

very memory inefficient when ?? is large, the method is restarted every so often and x?? reset to be x0

.

An incomplete GMRES algorithm function PerformGMRESRestarted() has been written in

linear_algebra.cpp.

A key component of the GMRES algorithm is the Arnoldi iteration that seeks to find an orthonormal basis

of ????. At the ??th step of the iteration, the Arnoldi method constructs the following matrix decomposition

of ??:

?????? = ????+1???

??,

where the columns of ???? (????+1) contain the orthonormal basis of ???? (????+1, resp.) and ???

?? is a (??+1)× ??

upper Hessenberg matrix. That is, a matrix that is nearly upper triangular but has non-zero components

on the first subdiagonal.

The ??th step of the Arnoldi algorithm is:

Algorithm 1 One step of the Arnoldi Iteration.

Require: ?? > 0, ??, ????:

1: Let q?? be the ??th column of ????.

2: Let h = {???

}

??+1

??=1 be a vector of length ?? + 1.

3: Compute q??+1 = ??q??

4: for ?? = 1, … , ?? do

5: ??? = q??+1 ? q??

.

6: q??+1 = q??+1 ? ???q??

.

7: end for

8: ???+1 = ‖q??+1‖2

.

9: q??+1 = q??+1/???+1.

10: ????+1 = [????, q??+1].

11: return ????+1 and h.

(a) In linear_algebra.cpp, write a C++ function which implements one step of the Arnoldi iteration

method defined above.

The function should have the following prototype

void PerformArnoldiIteration ( csr_matrix & matrix ,

dense_matrix & krylov_matrix , int k, double* hessenberg )

MATH4063 Turn Over

4 MATH4063

Here, matrix is ??, k is the step of the iteration to perform, krylov_matrix is the matrix containing

the orthonormal basis, where each row is a basis vector. Upon entry, krylov_matrix should have ??

rows and upon exit it should contain ?? + 1 rows, with the new basis vector in the last row.

Finally, upon exit, hessenberg should contain h, which is the final column of ???

??. You may assume that

hessenberg has been preallocated to be of length ??+1 before the call to PerformArnoldiIteration.

Your function should make use, where possible, of prewritten functions defined in dense_matrix.cpp

and vector.cpp. Your code should also make use of the matrix multiplication function from Q1.

Once you have written PerformArnoldiIteration() the GMRES function should function as intended.

Note: Storage of the basis functions in the rows of krylov_matrix, rather than in the columns,

improves efficiency of the code.

(b) In csr_matrix.cpp, write a C++ function that will read from a file a matrix already stored in CSR

format and a vector. You may assume the file structures are as in matrix1.dat and vector1.dat on

Moodle and you may use these data files to test your function.

(c) Write a test program in file q2c.cpp that will read in the matrix ?? from matrix2.dat and the vector

x from vector2.dat, compute b = ??x, then use PerformGMRESRestarted() with the default input

arguments to find an approximation x?to x. At the end of the calculation, print to the screen the error

‖x ? ? x‖2

.

[20 marks]

3. The file mesh.hpp contains a struct that defines a mesh data structure mesh for a general mesh comprising

axis-aligned rectangular cells. In particular, each cell in the mesh has an additional struct called

cell_information that contains, among other things, information about the cell neighbours. Familiarise

yourself with these data structures by looking in mesh.hpp.

mesh.cpp contains two functions that will generate meshes, they are:

? ConstructRectangularMesh() - this constructs a mesh on the rectangular domain ???? = [??, ??] ×

[??, ??].

? ConstructLShapedMesh() - this constructs a mesh on the L-shaped domain ???? = ????\????, where

???? = [(?? + ??)/2, ??] × [(?? + ??)/2, ??].

(a) In finite_volume.cpp, write a C++ function that will create the storage for a matrix ?? in CSR format

and a RHS vector F required for a cell-centred finite volume method for solving (1)-(2). You should

follow the procedure outlined in the Unit 6 lecture notes. As one of the inputs, your function should

take in a variable of type mesh.

(b) In csr_matrix.cpp, write a C++ function that will output to the screen a matrix stored in CSR format

in the same style as in matrix1.dat.

(c) In Q3c.cpp, write a program that will ask the user to supply the number of cells in each coordinate

direction of a rectangular mesh, sets up the mesh using ConstructRectangularMesh() then calls the

function from part (a) to set up the corresponding matrix and finally prints it to the screen using the

function from part (b).

[20 marks]

MATH4063

5 MATH4063

4. (a) In finite_volume.cpp, write a function that takes in a mesh, uses the function from Q3(a) to construct

?? and F, then populates it with the correct entries to solve problem (1)-(2) using the cell-centred finite

volume method, as outlined in the Unit 6 notes. The function should also take as input the functions

??(??, ??), b(??, ??) and the Dirichlet boundary function ??(??, ??).

(b) In Q4b.cpp, write a main program to ask the user to select from the following problems and supply

the number of cells in each coordinate direction.

1. ? Rectangular Mesh - ?? = 0, ?? = 1, ?? = 0 and ?? = 1;

? ??(??, ??) = 1;

? ??(??, ??) = 0;

? b = 0.

2. ? L-shaped Mesh - ?? = 0, ?? = 1, ?? = 0 and ?? = 1;

? ??(??, ??) = 8??2

cos(2????) cos(2????);

? ??(??, ??) = cos(2????) cos(2????);

? b = 0.

3. ? Rectangular Mesh - ?? = ?1, ?? = 1, ?? = ?1 and ?? = 1;

? ??(??, ??) = 1;

? ??(??, ??) = 0;

? b = (10, 10)?.

4. ? L-Shaped Mesh - ?? = 0, ?? = 1, ?? = 0 and ?? = 1;

? ??(??, ??) = 0;

?

??(??, ??) = {

1, ?? = 0, 0.25 < ?? < 0.75,

0, otherwise;

? b = (

50??

√??2+??2

,

?50??

√??2+??2

)

?

.

The code should then set up the linear system arising from the finite volume discretisation and solve

the system

??u? = F

using PerformGMRESRestarted().

Finally, print to the screen the maximum value of u?.

Hint: Once you have computed u? you can output it to together with the mesh to a file using

OutputSolution() in finite_volume.cpp. plot_solution.py can then be used to plot the solution in

Python.

Note, if you are unable to get the iterative solver from Q2 working, then you may create the finite volume

matrix ?? as if it were a dense matrix (i.e store all the zero entries) and use the function

PerformGaussianElimination() from dense_matrix.cpp to solve the system of equations. This will incur

a small penalty. Note, an illustration of the use of PerformGaussianElimination() can be found in the

main program inside gaussian_elimination_test.cpp.

[30 marks]

MATH4063 End