Updated on March 4, 2023 Spring 2022 Math 3607: Project 2 Due: 8:00PM, Friday, March 10, 2023

欢迎使用51辅导,51作业君孵化低价透明的学长辅导平台,服务保持优质,平均费用压低50%以上! 51fudao.top

Updated on March 4, 2023

Spring 2022 Math 3607: Project 2 Due: 8:00PM, Friday, March 10, 2023

 Instructions

• You are not allowed to use MATLAB commands and functions other than the ones discussed

in lectures, accompanying live scripts, textbooks, and homework/practice problem solutions.

• You may be requested to explain your code to me, in which case a proper and satisfactory explanation must be provided to receive any credits on relevant parts.

• You are allowed to collaborate in a group of up to three students. In such a case, submit a single file to Gradescope. The uploader must select all group members at the time of upload. If the uploader does not select the group members, they will not receive credit.

• If any code is found to be plagiarized from the internet or from other students, you will receive a zero on the entire project and will be reported to the COAM.

• Do not carry out computations using Symbolic Math Toolbox. Any work done using sym, syms, vpa, and such will receive NO credit.

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

• Answers to analytical questions (ones marked with à ) without supporting work or justifi- cation will not receive any credit.

                1

 

1 PLU Factorization via Outer Product Elimination [30 points]

In class, we derived a PLU factorization algorithm using outer product elimination. We wrote

n

ÿT lkuk ,

k“1

where lk is the kth column of L and ukT is the kth row of U. Recall that L is a row-permuted unit lower triangular matrix with i “ pi1, i2, . . . , inq encoding the row-permutation history. Below is the summary of the steps:

ˇˇ››

Step 1. Define A1 “ A. Find the index i1 such that ˇApi1, 1qˇ “ ›Ap:, 1q› 8. Then set

u 1T Ð A 1 p i k , : q

l1 Ð A1p:, 1q{u1,1

%PLUFACT

% L,U %

% i

% s

PLU factorization.

row-permuted unit lower triangular and upper triangular

matrices such that A = LU.

permuation vector such that L(i,:) is a unit lower triangular

number of row swaps

end end

n = length(A);

L = zeros(n);

U = zeros(n);

i = zeros(1,n);

s = 0;

for k = 1:n

    %% [PART 1] Determine i_k and update s.

    % (An if-statement is allowed.)

A “ LU “

Tˇˇ Step2. Fork“2,...,n,defineAk “Ak ́1 ́lk ́1u . Findtheindexik suchthatˇApik,kqˇ“

(a)

›› k ́1 ›Ap:, kq› 8. Then set

u kT Ð A k p i k , : q

lk Ð A1p:, kq{uk,k

Complete the following MATLAB function by implementing the algorithm above.

function [L,U,i,s] = plufact(A)

 %% [PART 2] Set U(k,:) and L(:,k), then update A.

% (A for-loop is NOT allowed.)

 2

 

Note 0. All you need to do is to fill in the two parts inside the main loop. Do not modify anything else. In the loop, you are allowed to use one if-statement, but no additional for-loop is allowed. This program carries out PLU factorization and counts the number of row swaps. Include the function at the end of your live script.

Note 1. This is different from myplu from the textbook or lecture slides.

Note 2. s keeps track of the number of row swaps, which means that s needs to be incremented

by 1 each time ik is different from k.

Note 3. Each part inside the main loop should only take about three to five lines of code. If

you are writing anything longer, then it is very likely that you are doing something wrong.

Note 4. Don’t waste your time googling. Review the code I wrote in class or lufact in “Operation Counts and PLU Factorization” live script.

Tip. You may use the following code block to test your code. Each of the three checks should yield zero.

A=[204 3 -2 0 2 -13 1 15 2 -4.5

-4 5 -7 -10];

[L,U,i,s] = plufact(A);

[L1,U1,i1] = lu(A,’vector’); % MATLAB’s PLU factorization routine L_check = norm(L(i,:)-L1)

U_check = norm(U-U1)

i_check = norm(i-i1)

(b) à Explain why, if P A “ LU is a PLU factorization, n

(c)

(d)

            uii,

Using the previous part, write a MATLAB function mydet2 that computes the determinant of a given matrix A using plufact function written for part (a). Include the function at the end of your live script.

Use your function and the built-in det on the matrices gallery(’cauchy’, n) for n “ 3, 4, . . . , 8, and make a table using fprintf showing n, the determinant calculated using your function mydet2, and the relative error when compared to det.

detpAq “ p ́1q

where s is the number of row swaps needed to transform I to P.

sź i“1

    3

 

2 Visualization of Matrix Norms in 3-D [30 points]

Recall that

}A}p “ max }Ax}p , p P r1,8s. }x}p “1

In this problem, we generate three-dimensional visualization of this definition.

(a) Complete the following program which, given p P r1,8s and A P R3ˆ3, approximates }A}p and plots the unit sphere in the p-norm and its image under A. Avoid using loops as much as possible. Include the function at the end of your live script.

    function norm_A = visMatrixNorms3D(A, p)

        %% Basic checks

        if size(A,1)~=3 || size(A,2)~=3

            error(’A must be a 3-by-3 matrix.’)

        elseif p < 1

            error(’p must be >= 1.’)

end

        %% Step 1: Initialization

        nr_th = 41; nr_ph = 31;

        th = linspace(0, 2*pi, nr_th);

        ph = linspace(0, pi, nr_ph);

        [T, P] = meshgrid(th, ph);

        x1 = cos(T).*sin(P);

        x2 = sin(T).*sin(P);

        x3 = cos(P);

        X = [x1(:), x2(:), x3(:)]’;

        %% Step 2: [FILL IN] Normalize columns of X into unit vectors

        %% Step 3: [FILL IN] Form Y = A*X; calculate norms of columns of Y

        %% Step 4: [FILL IN] Calculate p-norm of A (approximate)

        %% Step 5: [FILL IN] Generate surface plots

end

The following steps are carried out by the program. Step 1. Create 3-vectors

»fi

cos θi sin φj

xk “—–sinθisinφjffifl, for1ďiď41,1ďjď31 (1)

cos φj

using 41 evenly distributed θi in r0,2πs and 31 evenly distributed φj in r0,πs. Note the

use of meshgrid, which is useful for surface plots later.

Step 2. Normalize xk into a unit vector in p-norm by xk Ð xk{}xk}p (that is, replacing xk with xk{}xk}p).

  4

 

Step 3. For each k, let yk “ Axk. Calculate and store }yk}p.

Step 4. Approximate }A}p based on the norms }yk}p calculated in the previous step.

Step 5. Generate surface plots of the unit sphere in the p-norm and its image under A. Use surf function. Use subplot to put two graphs side by side.

(b) Run the program with p “ 1, 32 , 2, 4, all with the same matrix »fi

200

A “ —–0 cospπ{12q  ́ sinpπ{12qffifl , (2)

0 sinpπ{12q cospπ{12q by executing the following code block.

 A = [2 0 0;

     0 cos(pi/12) -sin(pi/12);

     0 sin(pi/12)  cos(pi/12)];

visMatrixNorms3D(A, 1);

visMatrixNorms3D(A, 3/2);

visMatrixNorms3D(A, 2);

visMatrixNorms3D(A, 4);

% Depending on how you write the code,

% you may need to use ’clf’ or ’hold off’

% in between function calls.

      Figure 1: Example output.

5

 

 


51作业君

Email:51zuoyejun

@gmail.com

添加客服微信: abby12468