G12ISC代做、代写Scientific Computation、代做Matlab编程、Matlab代写

日期：2018-12-12 09:34

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.