关于程序员:解析-G12ISC-Scientific-Computation

G12ISC Introduction to Scientific Computation
— Coursework 2 (10%) —
Submission deadline: 4pm, Wednesday, 12 Dec 2018
This coursework contributes 10% towards the overall grade for the module.
Note: This Coursework PDF document is complete. A total of 40 points can be obtained
in this coursework. The weighting is enclosed in brackets for each question, e.g. [4 / 40].
See Moodle: Assessed Coursework 2 for all grading criteria per question.
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.
Coursework Aim: In this coursework you will develop Matlab code related to algorithms that
solve linear systems of equations.
Getting Started:
Download the contents of the“Coursework 2 Pack”folder from Moodle.
Do the coursework by simply adding MATLAB code to the downloaded file cw2.m, and by
completing the incomplete function m-files: backSub.m, forwElimPP.m, forwElimPPperm.m,
genIter.m, genIterTol.m.
Your responses to the numbered questions below should be written in cw2.m, except the
responses to questions 1, 4, 6, 8, 11 which involve the completion of the other m-files.
Advice:
Write your code as neatly and concisely as possible so that it is easy to follow.
Add some comments to your scripts that indicate what a piece of code does.
Remember to format output using for example disp, num2str or fprintf.
Some marking criteria: Correctness of numerical output, clarity of code and comments, completeness
of figures (title, axis labels, legend, etc) and tables, clear and concise answers to
open questions, . . .
Submission Procedure: The coursework is to be submitted using the“Assessed Coursework
2”assignment on Moodle. A complete submission consists of the following files:
All m-files (including those that are already complete): cw2.m, forwElim.m, backSub.m,
testMat.m, forwElimPP.m, forwElimPPperm.m, genIter.m, genIterTol.m
1 PDF file: cw2.pdf This PDF file is created using Matlab“Publish”on cw2.m (see for example
Section 1.2.9 of the Matlab Guide (by Tom Wicks) for instructions on publishing scripts).
I Gaussian elimination without pivoting
To solve a linear system of equations, consider the Gaussian elimination algorithm
consisting of forward elimination without row interchanges and backward substitution.
The function m-file forwElim.m already contains a complete implementation of
School of Mathematical Sciences — Version: 8th December, 2018 — Lecturer: Dr Kris van der Zee Page 1 of 5
the forward-elimination step (it will be helpful to familiarise yourself with its contents).
This m-file essentially implements Algorithm 6.1, steps 1, 4, 5 & 6, as discussed in
the Lectures (see Lecture 4 and 5, slides & notes, available on Moodle).
1 Modify the function m-file backSub.m. This function m-file should implement the [4 / 40]
backward-substitution step (which is essentially steps 8 & 9 of Algorithm 6.1),
and has the following prototype:
function x = backSub (A)
where A is an n × (n+1) augmented matrix in echelon form, and the output x is
the corresponding solution vector (which is a column vector with n entries).
2 By using the function m-file forwElim.m and your backSub.m, compute the ech- [2 / 40]
elon form and the solution x in case the system of linear equations has an augmented
matrix A1 given by
A1 =
1 1 2 1 8
1 1 1 0 2
2 2 3 3 20
1 1 4 3 4
Make sure your cw2.m outputs the augmented matrix in echelon form and the
solution vector x, and that this is visible when“publishing”cw2.m as a PDF.
Hint: You should verify that the solution is (column vector) x = [?7; 3; 2; 2].
3 Consider the matrix A of size n × n and vector b of size n × 1 as provided by the [4 / 40]
function [A, b] = testMat(n). Take a large range of values for n (for example,
n = 50, 100, 150, 200, . . . , 950, 1000, . . . up to some large number that depends
on how good your computer is, and how much patience you have), and use
Matlab to compute the time it takes for forwElim.m to obtain the echelon form
for the augmented matrix [A, b], as well as the time it takes for your backSub.m
to subsequently obtain the solution vector.
Plot the computed time for forwElim.m versus n in a figure using loglog, and
plot in the same figure the computed time for backSub.m.
Hint: Use tic and toc to set Matlab’s stopwatch.
Based on your plots, provide an explanation for the observed behaviour of the
timings.
I Partial pivoting Now consider Gaussian elimination with forward elimination using
partial pivoting.
4 Complete the function m-file forwElimPP.m. The function has the following pro- [4 / 40]
totype:
function B = forwElimPP (A)
where A is an augmented matrix of size n × (n+1), and the output B is the
corresponding echelon form obtained by doing forward elimination and partial
Coursework 2 Page 2 of 5
pivoting (see Lecture 5).
Hint: To find a maximum in some vector, you could use a for loop and boolean
comparisons, or you can use the Matlab function max (see help max for more
info).
5 Repeat Question 2, but now using your forwElimPP.m and in case of the aug- [2 / 40]
mented matrix A2 given by
A2 =
0 1 2 1 1
1 1 1 0 8
2 2 3 3 20
1 1 4 3 4
(Note: The solution is the same as before.)
6 Complete the function forwElimPPperm.m, which has the following prototype: [3 / 40]
function [B , P] = forwElimPPperm (A)
This function implements the same algorithm as forwElimPP.m, but additionally,
it outputs P, which is the n × n permutation matrix corresponding to the pivoting
process. In other words, P is a rearrangement of the n × n identity matrix, so
that PA is the row-rearrangement of A. Note that the forward elimination without
pivoting of PA gives exactly the same echelon form as forward elimination with
pivoting of A.
7 Apply your forwElimPPperm.m to A3, which is given in Matlab notation by [2 / 40]
A3 = […
26 -34 31 -48 25 -18 13 -32 -28 13 29
0 0 25 42 -27 31 -14 23 -13 -48 67
24 39 -38 15 -44 29 0 -13 -41 41 -152
-39 2 3 43 27 0 -28 34 14 0 362
18 20 -17 -34 17 1 15 23 -32 25 292
-4 -35 5 42 22 14 10 7 -45 31 281
0 45 -10 29 14 45 -11 -32 22 -12 421
-40 4 -8 8 -8 -6 -36 46 -15 12 -65
32 18 -32 -6 -11 -44 0 0 0 0 -337
-32 -46 -24 -24 32 37 -8 42 0 3 124
];
(You should be able to select the above matrix and copy–paste it into Matlab.)
(Alternatively, find an m-file with this matrix on the discussion forum on Moodle.)
Output the permutation matrix P.
Apply forwElim.m to A = PA3, and output the resulting echelon form.
I Jacobi method and Gauss–Seidel method
To solve a linear system of equations, consider iterative methods that can be written
Coursework 2 Page 3 of 5
as the following general iteration method:
x
(k) = Tx
(k1) + c for k = 1, 2, . . . , Nmax.
8 Complete the function m-file genIter.m, which has the following prototype: [4 / 40]
function x_mat = genIter (T ,c , x0 , Nmax)
where the input consists of an (n × n) iteration matrix T, an (n × 1) vector c,
the (n × 1) vector x0 (a given initial approximation), and Nmax, which is the
total number of iterations to be peformed by the method. The output x mat is
a matrix of size n × (Nmax+1) that stores all approximation vectors xk for k =
0, 1, . . . , Nmax, that are generated by the iteration method.
9 Consider the linear system Ax = b, where [3 / 40]
A =
4 1 0
1 8 1
0 1 4
and b =
48
12
24
.
By using your genIter.m, obtain the approximations as provided by the Jacobi
method and by the Gauss–Seidel method, taking in both methods x0 = [1; 1; 1]
and Nmax = 12. (The output can be suppressed.)
Hint: You can simply pre-compute the T and c as used in these methods (see
Lecture 7, slides & notes, available on Moodle).
Create (and output) for each method a table with a column for k (ranging from
0 to Nmax), and three columns for the corresponding approximations x
(k)
1
, x
(k)
2
,
and x
(k)
3
. (Hence, two tables with four columns each are requested.)
10 The exact solution is x = [13; 4; 7]. [4 / 40]
Create (and output) a semilogy figure with a plot of the l2 norm of the error,
kx x
(k)k2, versus k for each method. Also do this for the l∞ norm of the error,
kx x
(k)k∞, for each method. (Hence, two figures with two plots (one for Jacobi
and one for Gauss–Seidel) are requested.)
Hint: Type“help norm”in Matlab to see how Matlab can compute the l2
norm k · k2 and the l∞ norm k · k∞ of a vector.
From your semilogy plots, what is (approximately) the observed slope of the
results? (Four values are requested, one for each plot.)
Following questions available since 4 Dec 2018. . . (= last set of questions)
I Jacobi method and Gauss–Seidel method with stopping criterion
11 To improve on genIter.m, complete the function m-file genIterTol.m, which [2 / 40]
has the following prototype:
function [x_mat , Niter , resNorm_vec] = genIterTol (T ,c , x0 , Nmax ,
tol ,A , b )
Coursework 2 Page 4 of 5
where the input consists additionally of a (positive) real number tol, the n × n
matrix A and n × 1 vector b. This m-file implements the general iteration method
with a stopping criterion based on the l2 norm of the residual vector, that is, the
iterations are stopped if
b Ax
(k)2
< tol . The output consists additionally
of Niter, which is the total number of iterations required to reach the stopping
criterion (or to reach Nmax, whichever comes first), and resNorm vec which is a
vector of size Niter and stores all values of
b Ax
(k）2
, k = 1, 2, . . . , Niter.
Note that the output x mat is now of size n × (Niter+1).
Hint: Use the command“break”to prematurely exit“for”or“while”loops.
12 Consider again the linear system defined in Question 9. By using [2 / 40]
your genIterTol.m, obtain the approximations as provided by the Jacobi
method and by the Gauss–Seidel method, taking in both methods x0 = [1; 1; 1]
and tol = 1e?10. Choose Nmax sufficiently large, so that Niter is determined
by the stopping criterion. (The output can be suppressed.)
Create (and output) for each method a table with a column for k (ranging from 1
to Niter) and a column for resNorm vec. (Hence, two tables with two columns
each are requested.)
13 Consider the Jacobi method for the linear system Ax = b, where A and b are [4 / 40]
the n × n matrix and n × 1 vector defined by:
A =
2+δ 1 0 · · · 0
1 2+δ 1
0 · · · 0 1 2+δ
and b =
, respectively.
where δ ≥ 0. (Matrix A appears in the numerical solution of certain differential
equations.)
Use Matlab to output the matrix T and vector c for the Jacobi method when
δ = 0.5 and n = 7.
Assuming δ = 0, use your genIterTol.m to study how the total number of iterations
Niter for the Jacobi method depends on n as n → ∞. Take tol = 1e3
and x0 = [0 ; . . . ; 0] (and Nmax sufficiently large).
Finally, study how Niter depends on δ, for fixed large values of n (e.g.,
n = 50). What happens as n → ∞ Take tol and x0 as before. It is suffi-
cient to consider values of δ that are ≥ 0.
Hint: Provide sufficient evidence (for example output and/or figures) to support
your answers. You are allowed to create additional function m-files.