Math551 - lab06 - Solved

Goal: In this assignment you will learn basic Matlab commands dealing with LU factorization, learn to find LU factorizations of matrices, and learn to find inverses and solve linear systems using the LU factorization. During the lab session, your lab instructor will teach you the necessary MATLAB commands to complete the assignment.


Variables: x1, A, q
At times you will also be asked to answer questions in a comment in your M-file. You must format your text answers correctly. Questions that you must answer are shown in gray boxes in the lab. For example, if you see the instruction

Q7: What does the Matlab command lu do? your file must contain a line similar to the following somewhere

% Q7: It computes the LU decomposition of a matrix.

The important points about the formatting are

•    The answer is on a single line.

•    The first characters on the line is the comment character %.

•    The answer is labeled in the form Qn:, where n stands for the question number from the gray box.

If you do not format your answers correctly, the grading software will not find them and you will lose points.

Tasks
Create an M-file named m551lab06.m in which you will complete the following tasks. 1. Start your M-file by defining the variable A to hold the matrix

A   .

Variables: A
2. Add the command

[LA,UA]=lu(A);

Variables: LA, UA
Execute the M-file and examine the contents of the matrices LA and UA. Observe that UA is upper triangular.

Q1: Is LA lower triangular? Why?

Hint: Type doc lu and look at the difference between the commands

                                             [L,U] = lu(A)                                     and                                       [L,U,P] = lu(A)

3. Now add the command

[L,U,P]=lu(A);

Variables: L, U, P
Execute the M-file again and look at the matrices L, U and P.

Q2: Is L lower triangular?

4.    In the Matlab command window, type the commands

>> P’*P

>> P*P’

Q3: What do the results tell you about the relationship between P and PT?

5.    In the Matlab command window, type



 
 
 

>> A-P’*L*U

and examine the result. It should be consistent with what you saw in the lu documentation. (Round-off errors may prevent you from seeing exactly what you expect, but it should be close.)

6.    From the previous step we saw that A = PTLU. For your own edification, you should check that this means that A−1 = U−1L−1P. One way to compute the matrix on the right-hand side is to use the backslash operator: U\(L\P). Add a variable named inv  LU to your M-file, storing the inverse of A computed in this way.

Variables: inv LU

7.    Add a command to compute the inverse by applying the inv function to A. Name the result inv A.

Variables: inv A
8.       Define a variable called inv err holding the difference between inv  LU and inv A. Variables: inv err

Q4: Is inv err the zero matrix? Why not?

9.       Create a variable b holding the the column vector

b   .

Variables: b

10.   Now we’re going to solve the linear system

Ax = b
using the LU decomposition. This can be done efficiently in Matlab by observing that

                                                Ax = b         ⇐⇒          PTLUx = b          ⇐⇒        x  .

We’ll work our way from the innermost parentheses outward and use the backslash operator instead of explicitly computing inverses. As we proceed, notice that we never use the expensive operations of inv or matrix-matrix multiplication.

Now, add the line

Pb = P*b; beneath the definition of b in your M-file. Variables: Pb

Run the M-file and compare b and Pb in command window.

Q5: How are the two vectors related?



 
 
 

11. Add two more lines to your M-file, one computing L\Pb and storing it in the variable y, and the next finishing the computation by computing U\y and storing it in x.

Variables: x, y
Run your M-file, and examine the variable x. If you have completed all steps correctly, x should be the solution to the linear system.

>> x x =

1.0000

2.0000

3.0000

If the above isn’t what you see, go back and check your work.

12. Now define a new variable z in your M-file and assign it the value A\b and define xz err to be the difference between x from the previous steps and z.

Variables: z, xz err
In the command window, examine the value of xz err. This time you should get exactly the zero vector. The reason is that the steps you took in finding x are exactly the same steps that Matlab took in finding z.