MTH2051 Introduction to Computational Mathematics Assignment 2. Numerical linear algebra

SCHOOL OF MATHEMATICAL SCIENCES

MTH2051/3051

Introduction to Computational Mathematics

Assignment 2. Numerical linear algebra.

Due date: Monday August 27th, 3pm. (Week 4.)

Late submissions will be subject to late penalties (see the unit guide for full details).

Please submit the hardcopy of your completed assignment in the assignment box of your tutorial class. A printout of the Matlab code for the programming exercises as well as relevant output are to be submitted as part of your hardcopy submission. In addition, collect all your Matlab m-files in one zip or rar-file and submit this file through Moodle. The Moodle submission is due at the same time as the hardcopy submission.

1. LU decomposition. (10 marks)

Find the LU decomposition of

A =

 1 4 72 5 8 3 6 10

 . (1) Show all your work, including explicit expressions for the matrices Mi and Li.

2. Computational cost of computing a matrix determi-

nant. (20 marks)

We know from class that computing the determinant of an n × n matrix using LU decomposition has computational cost W = O(n3). But how expensive is it to compute a determinant using the common recursive definition formula for determinants?

(a) Determine the number of floating point operations required for calculating the determinant of an n × n matrix using a straightforward implementation of the recursive definition of the determinant. (You can combine additions and multipli- cations.) Hint: consider the number of flops done on each recursive level separately. How

School of Mathematical Sciences Monash University

many determinants are to be computed on each recursive level? On each recursive level, how many flops does it take to compute each of these determinants? (By us- ing determinants that are one size smaller, computed on the next recursive level.) Then sum up over all the levels (give the result as a sum like

∑n i=2 . . .).

(b) Find an approximation for the expression derived that is valid for large n. (Hint: this approximate result should be a very simple expression, proportional to n!, and you will need expx = 1 + x + x2/2! + x3/3! + x4/4! + . . .).

3. LU decomposition with partial pivoting. (20 marks)

Solving a linear system A~x = ~b by LU decomposition requires three phases: compute the decomposition A = LU , solve L~y = ~b (forward substitution), and solve U~x = ~y (backward substitution).

Download the m-file LUFactorization.m from Moodle. This program implements the LU decomposition pseudocode from class, not worrying about pivoting (it checks for a zero pivot element and stops). The program returns the matrices L and U as one matrix Q, using the compressed storage method discussed in class: L is stored in the part of Q below the diagonal (the 1s on the diagonal of L are not stored in Q), and U is stored in the upper triangular part of Q.

(a) You are asked to complete the implementation and code up forward and backward substitution in Matlab, using only primitive commands (i.e. do not appeal to the built-in linear algebra functionality in MATLAB). You should use the function header given for each below and submit a copy of your code on Moodle and in the assignment package submitted to the assignment boxes.

(i) Forward Substitution. This function should take as input the matrix Q output by LUFactorization.m.

function y = ForwardSubstitution(Q, b)

% Usage: y = ForwardSubstitution(Q, b)

% Perform forward substitution using the unit

% lower triangular portion of the matrix Q.

(ii) Backward Substitution. This function should take as input the matrix Q output by LUFactorization.m.

function x = BackwardSubstitution(Q, y)

% Usage: x = BackwardSubstitution(Q, y)

% Perform backward substitution using the

% upper triangular portion of the matrix Q.

(b) Download the m-file VerifyLU.m from Moodle to test your implementation using randomly generated matrices. Choose n = 123, run VerifyLU(n) and provide n and the reported error you obtain in your assignment package. Clearly, the error you report must be the error you obtained when running the code you submit in your assignment package.

2018 2

School of Mathematical Sciences Monash University

(c) The LU decomposition code you downloaded cannot handle zero pivot elements, so it cannot solve simple linear systems like[

0 1 1 1

] ~x =

[ 1 2

] . (2)

To fix this, you are asked to modify LUFactorization.m and implement LU de- composition with partial pivoting. In this variant of LU decomposition, you choose the pivot element in matrix A(k) to be maxi≥k |a(k)ik |, i.e., you choose as pivot the matrix element in column k at or below a

(k) kk that is the largest in absolute value.

This means that, in every step of the LU decomposition, you may have to swap two rows in the U and L you have computed so far, and you have to keep track of these swaps. It is easiest to keep track of the swaps in a permutation matrix P (which can also be obtained by swapping rows, starting from the identity matrix

I). After computing PA = LU , you then solve L~y = P~b (forward substitution), followed by U~x = ~y (backward substitution).

This PA = LU decomposition with partial pivoting handles the case of zero pivot elements that may occur in the usual LU decomposition, and it turns out that always choosing the largest pivot element in column k also tends to lead to more accurate numerical results than the usual LU decomposition. (We will discuss in class that this variant is more numerically stable.)

You should use the function header given below for your PA = LU decomposi- tion and submit a copy of your code on Moodle and in the assignment package submitted to the assignment boxes.

function [P, Q] = PALUFactorization(A)

% Usage: [P, Q] = PALUFactorization(A)

% Compute the PALU factorization of the matrix A

% and return P, and L and U stored in the single matrix Q.

(d) Download the m-file VerifyPALU.m from Moodle to test your PALU implementa- tion. How much accuracy gain do you see when using partial pivoting versus no pivoting, for a random matrix with n = 201?

Question 5 below is only for students enrolled in MTH3051 (not for MTH2051 students).

2018 3

School of Mathematical Sciences Monash University

4. Vandermonde matrices and ill-conditioning. (10 marks)

(Only for students enrolled in MTH3051; not for MTH2051

students.)

A Vandermonde matrix is an n×n matrix formed from a vector ~w = (w0, w1, w2, . . . , wn−1) as follows:

V (~w) =

 wn−10 wn−11 wn−12

… wn−1n−1

wn−20 wn−21 wn−22

… wn−2n−1

· · ·

w20 w21 w22 …

w2n−1

w0 w1 w2 …

wn−1

1 1 1 … 1

 . (3) You may generate Vandermonde matrices using the MATLAB command vander, which takes the vector ~w as its single parameter.

(a) Consider the Vandermonde matrix generated with ~w = [10.0 : 0.1 : 11.0]. Solve

the system V (~w)~x = ~b where ~b = [1; 9; 1; 9; 1; 9; 1; 9; 1; 9; 1]T using LU decomposi- tion and Forward/Backward solver. Plug your answer back into the system and compute V (~w)~x. You would expect to get a result that is close to the original

right-hand side ~b, but what do you observe? Provide the result you computed. (You can use your code from Question 3, or Matlab’s built-in linear solver. You can check (not required) that the partial-pivoting variant PALU does in this case not lead to more accurate results than basic LU .)

(b) Compute the condition number of the matrix V (~w) and use this to explain the observed discrepancy. (You may calculate the 2-norm condition number of a ma- trix A using the cond command. This command returns an explicit value of the condition number for the provided matrix.) Note: You will see later in this unit that solving a Vandermonde system can in prin- ciple be used to compute the coefficients of an interpolating polynomial, but this example clearly illustrates that solving a Vandermonde system using LU decom- position is often a bad idea due to the ill-conditioning of Vandermonde matrices. (Lagrange or Newton polynomials are a much better approach.)

Submit a copy of your Matlab code for items (a) and (b) on Moodle and in the assignment package submitted to the assignment boxes.

5. What to hand in

When submitting your hardcopy assignment package to the assignment box of your tutorial class you should include

F Complete written or typed answers to the questions.

2018 4

School of Mathematical Sciences Monash University

F Results and plots generated by your Matlab programs, and printouts of your Mat- lab code.

In addition, make sure to also submit your Matlab programs on Moodle. Marks will be awarded for code correctness (also make sure to include well-chosen comments in your code; it should be easy for anybody to read and understand your program; also, please use proper indentation in your code).

2018 5