Begin Sci Comp5

Updated on September 28, 2021

Math 3607: Homework 5 Due: 10:00PM, Tuesday, October 5, 2021

TOTAL: 30 points

• Problems marked with v are to be done by hand; those marked with � are to be solved using a computer.

• Important note. Do not use Symbolic Math Toolbox. Any work done using sym or syms will receive NO credit.

• Another important note. When asked write a MATLAB function, write one at the end of your live script.

1. (Improved triangular substitutions; adapted from FNC 2.3.5) � If B ∈ Rn×p has columns b1, . . . , bp, then we can pose p linear systems at once by writing AX = B, where X ∈ Rn×p whose jth column xj solves Axj = bj for j = 1, . . . ,p:

A

 x1 x2 · · · xp

 

︸ ︷︷ ︸ =X

=

 b1 b2 · · · bp

 

︸ ︷︷ ︸ =B

.

(a) Modify backsub.m and forelim.m from lecture1 so that they solve the case where the second input is an n×p matrix, for p ≥ 1. Include the programs at the end of your live script.

(b) If AX = I, then X = A−1. Use this fact to write a MATLAB function ltinverse that uses your modified forelim to compute the inverse of a lower triangular matrix. Include the program at the end of your live script. Then test your function using the following matrices, that is, compare your numerical solutions against the given exact solutions.

L1 =

 2 0 08 −7 0

4 9 −27

  , L−11 =

 

1 2 0 0 4 7 −

1 7 0

50 189 −

1 21 −

1 27

 

L2 =

 

1 0 0 0 1 3 1 0 0 0 13 1 0 0 0 13 1

  , L−12 =

 

1 0 0 0 −13 1 0 0

1 9 −

1 3 1 0

− 127 1 9 −

1 3 1

 

1Lecture 12 (09/24/21, F): Square Linear Systems (Introduction)

1

2. (Triangular substitution and stability; FNC 2.3.6) Consider the following linear system Ax = b: 

 1 −1 0 α−β β 0 1 −1 0 0 0 0 1 −1 0 0 0 0 1 −1 0 0 0 0 1

 

︸ ︷︷ ︸ A

  x1 x2 x3 x4 x5

 

︸ ︷︷ ︸ x

=

  α 0 0 0 1

  .

︸ ︷︷ ︸ b

(a) v Show that x = (1, 1, 1, 1, 1)T is the solution for any α and β. (b) � Using MATLAB, solve the system with α = 0.1 and β = 10, 100, . . . , 1012, making a

table of the values of β and |x1 − 1|. Write down your observation.

3. (Vectorizing mylu.m; FNC 2.4.7) Below is an instructional version of LU factorization code presented in lecture. function [L,U] = mylu(A) % MYLU LU factorization (demo only--not stable!). % Input: % A square matrix % Output: % L,U unit lower triangular and upper triangular such that LU=A

n = length(A); L = eye(n); % ones on diagonal % Gaussian elimination for j = 1:n-1 for i = j+1:n

L(i,j) = A(i,j) / A(j,j); % row multiplier A(i,j:n) = A(i,j:n) - L(i,j)*A(j,j:n);

end end U = triu(A);

end

Consider the innermost loop. Since the different iterations in i are all independent, it is possible to vectorize this group of operations, that is, rewrite it without a loop. In fact, the necessary changes are to delete the keyword for in the inner loop, and delete the following end line. (You should also put a semicolon at the end of i = j+1:n to suppress extra output.)

(a) � Make the changes as directed and verify that the function works properly. (b) v Write out symbolically (i.e., using ordinary elementwise vector and matrix notation)

what the new version of the function does in the case n = 5 for the iteration with j = 3.

4. (Application of LU factorization: FNC 2.4.6) When computing the determinant of a matrix by hand, it is common to use cofactor expansion and apply the definition recursively. But this is terribly inefficient as a function of the matrix size.

(a) v Explain why, if A = LU is an LU factorization,

det(A) = u11u22 · · ·unn = n∏

i=1 uii.

2

This part is an analytical question. Do it by hand. (b) �Using the result of part (a), write a MATLAB function determinant that computes

the determinant of a given matrix A using mylu from lecture. Include the function at the end of your live script. Use your function and the built-in det on the matrices magic(n) for n = 3, 4, . . . , 7, and make a table (using disp or fprintf) showing n, the value from your function, and the relative error when compared to det.

5. (Proper usage of lu; FNC 2.6.1) v Suppose that A ∈ Rn×n and b ∈ Rn. On the left is correct MATLAB code to solve Ax = b; on the right is similar but incorrect code. Explain using mathematical notation exactly what vector is found by the right-hand version. [L,U] = lu(A); x = U \ ( L \ b );

[L,U] = lu(A); x = U \ L \ b;

6. (FLOP Counting) v Do LM 10.1–12(a,b,d). Justify your calculation of p and c for each part.

7. (Matrix norms; Sp20 midterm) Let

A = [ 1 2 0 3

] .

(a) v Calculate ‖A‖1, ‖A‖2, ‖A‖∞, and ‖A‖F all by hand. (b) � Imagine that MATLAB does not offer norm function and you are writing one for

others to use, which begins with

function MatrixNorm(A, j) % MatrixNorm computes matrix norms % Usage: % mat_norm(A, 1) returns the 1-norm of A % mat_norm(A, 2) is the same as mat_norm(A) % mat_norm(A, ’inf’) returns the infinity-norm of A % mat_norm(A, ’fro’) returns the Frobenius norm of A

Complete the program. (Hint: To handle the second input argument properly which can be a number or a character, use ischaracter and/or strcmp.)

3