MAT 343 Arizona State University LU Factorization Laboratory 3 Exam 1. Differential equation exam will start in 25.5 hours and you will have 12 hours to fi

MAT 343 Arizona State University LU Factorization Laboratory 3 Exam 1. Differential equation exam will start in 25.5 hours and you will have 12 hours to finish it2. Applied linear algebra quiz will start in 22 hours and you will have 20 hours to finish it3. Differential equation Matlab deadline is in 60 hours4. Applied linear algebra Matlab deadline is in 122 hours MAT 343 Laboratory 3
The LU Factorization
In this laboratory session we will learn how to
1. Find the LU factorization of a matrix using elementary matrices
2. Use the MATLAB command lu to find the LU factorization
3. Solve linear systems using the LU factorization
Instructions: For your lab write-up follow the instructions of LAB 1.
Elementary Matrices and Row Reduction
An elementary matrix is any matrix that results from applying a single elementary row operation to the identity matrix. In this section we will discover some properties of elementary
matrices and their connection to row reduction of matrices.
For example, the following three elementary matrices were obtained from the 4 × 4 identity
matrix:

 
 

1 0 0 0
−5 0 0 0
0 0 1 0
0 1 0 0  0 1 0 0 −3 1 0 0

 
 

(1)
1 0 0 0  0 0 1 0  0 0 1 0
0 0 0 1
0 0 0 1
0 0 0 1
(The row operations were, respectively: interchange rows 1 and 3; multiply row 1 by −5; multiply row 1 by -3 and add to row 2).
We can construct each of the three types of elementary matrices using basic MATLAB
commands.
Elementary Matrices of Type I:
To swap rows i and j of the n × n identity matrix, set
E1 = eye(n);
E1([i,j],:) = E1([j,i],:)
Elementary Matrices of Type II:
To multiply row i of the n × n identity matrix by the scalar c, set
E2 = eye(n);
E2(i,i) = c
Elementary Matrices of Type III:
To replace row j of the n × n identity matrix by the sum of row j plus c times row i, set
E3 = eye(n);
E3(j,i) = c
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
1
Example 1
We can generate the three matrices listed in (1) using the following MATLAB commands.
>> E1 = eye (4);
>> E1 ([1 ,3] ,:)= E1 ([3 ,1] ,:)
E1 =
0
0
1
0
0
1
0
0
1
0
0
0
0
0
0
1
>> E2 = eye (4);
>> E2 (1 ,1)= -5
% generate the 4 x4 identity matrix
% swap rows 1 and 3
% generate the 4 x4 identity matrix
% change the entry in position (1 ,1) to -5
% equivalent to multiply row 1 by -5
E2 =
-5
0
0
0
0
1
0
0
>> E3 = eye (4); %
>> E3 (2 ,1)= -3 %
%
E3 =
1
0
-3
1
0
0
0
0
0
0
1
0
0
0
0
1
generate the 4 x4 identity matrix
change the entry in position (2 ,1) to -3
equivalent to multiply row 1 by -3 and add to row 2
0
0
1
0
0
0
0
1
Important Property: Suppose that E is an n × n elementary matrix obtained from the
identity by performing one of the elementary row operations. If A is an n × r matrix, then the
matrix EA is obtained by performing the same row operation on A.
EXERCISE 1
If you haven’t already done so, enter the commands in Example 1 to generate the elementary
matrices E1, E2 and E3 (you can suppress these matrices in your lab report).
Generate a random 4 × 3 matrix A (with integer entries), by typing
A = floor(10*rand(4,3))
Compute the products E1*A, E2*A and E3*A and compare A with each of the products. Describe
the effect of each multiplication on the matrix A. Be specific. For each of the three products,
describe specifically how the rows of A have changed.
What general pattern do you see? That is, what effect does multiplying a matrix A on the left
by an elementary matrix have on the matrix A?
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
2
The LU factorization using elementary matrices
EXERCISE 2
Enter the following matrix in MATLAB.


−2 1 −4
5 13 
A= 8
6 −9 17
(a) Determine elementary matrices E1 , E2 , E3 of Type III such that E3 E2 E1 A = U with U
an upper triangular matrix. The matrix E1 should turn the element in position (2,1) into
a 0. Enter this matrix in MATLAB as E1 using commands similar to the ones in Example
1. The matrix E2 should turn the element in position (3,1) into a zero. Enter this matrix
in MATLAB as E2. Note that to zero out the entries in column 1, you need to add or
subtract a multiple of row 1. Once you have found the matrices E1 and E2 , compute
the product E2 E1 A in MATLAB. Use format rat so that the entries will be given as
fractions. Based on the result, determine the matrix E3 that turns the element in position
(3,2) into a zero. Enter this matrix as E3 in MATLAB and compute U=E3*E2*E1*A.
(b) Compute the product L = E1−1 E2−1 E3−1 . The matrix L is lower triangular with ones on
the diagonal. Enter format short and verify that A = LU by computing A − LU in
MATLAB.
NOTE: From part (a) we have E3 E2 E1 A = U and therefore A = (E3 E2 E1 )−1 U = (E1−1 E2−1 E3−1 )U .
Since the elementary matrices and their inverses are lower triangular and have 1’s along the
diagonal, as is always true for elementary matrices of the third type, it follows that L is also
lower triangular with 1’s along the diagonal.
Permutation matrices
A permutation matrix is a square matrix that has exactly one 1 in every row and column and
0’s elsewhere. In this section we will look at properties of permutation matrices.
One way to construct permutation matrices is to permute the rows (or columns) of the
identity matrix. For example, we can construct


0 0 0 0 1
0 0 1 0 0



(2)
E=
0 0 0 1 0
1 0 0 0 0
0 1 0 0 0
by using the MATLAB commands
p = [5 ,3 ,4 ,1 ,2]; % permutation vector that defines
% the new order of the rows
E = eye ( length ( p ) ) ;
% define E as the identity matrix
E = E (p ,:) % permute the rows of E according to the
% permutation vector p
The second command creates the 5×5 identity matrix, and the third command uses the vector p
to permute its rows, so row p(1)= 5 becomes row 1, row p(2)= 3 becomes row 2, row p(3)= 4
becomes row 3 and so on (compare these row permutations with the vector p defined above).
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
3
EXERCISE 3
If you haven’t already done so, enter the commands in the example above to generate the
permutation matrix E defined in (2) (you can suppress this matrix).
Generate a 5 × 5 matrix A with integer entries using the command
A = floor(10*rand(5))
(a) Compute the product EA and compare the answer with the matrix A. How are the
two matrices related? Describe the effect on A of left multiplication by the permutation
matrix E. Be specific!
Compute the product AE and compare the answer with the matrix A. How are the two
matrices related? Describe the effect on A of right multiplication by the permutation
matrix E. Be specific!
(b) Compute E −1 and E T (recall that E T is computed in MATLAB with the command E’),
and observe that they are also permutation matrices. What else do you observe about
E −1 and E T ?
The LU factorization in MATLAB
Not all matrices have an LU factorization because not all matrices can be reduced to upper
triangular form using only row operation of Type III. However, if we allow the rows to be
interchanged, then all matrices have an LU factorization (but we have to keep track of the rows
we permute).
Example 2
The matrix


2 4 6
A= 1 2 5 
−1 5 3
does not have an LU factorization. Zeroing out the entries below the pivot in the first column
yields


2 4 6
 0 0 2 
0 7 6
and it is impossible to put this matrix in triangular form without permuting the rows. However,
the matrix


2 4 6
P A =  −1 5 3 
1 2 5
(obtained by interchanging row 2 and row 3 of
the permutation matrix

1

P = 0
0
A) does have an LU factorization. Here P is

0 0
0 1
1 0
In general, for any matrix A, it is always possible to find the LU factorization of P A, where
P is an appropriate permutation matrix.
The MATLAB command
[L,U,P] = lu(A)
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
4
returns a lower triangular matrix L (with 1’s on the diagonal), an upper triangular matrix U ,
and a permutation matrix P so that P*A = L*U.
For this example we have:
>> A=[2,4,6;1,2,5;-1,5,3]
A =
2
4
1
2
-1
5
>> [L,U,P]=lu(A)
L =
1
0
-1/2
1
1/2
0
U =
2
4
0
7
0
0
P =
1
0
0
0
0
1
6
5
3
0
0
1
6
6
2
0
1
0
and we can easily verify that P A = LU .
Remark: Note that, even if the LU factorization of A does exist, MATLAB will most likely
use a permutation matrix anyway. This is because MATLAB permutes the rows so that the
pivot is always the largest entry in the column. This technique, called partial pivoting, reduces
round-off errors.
Solving Systems using the LU factorization
Given A = LU , we can solve Ax = b in two steps:
(a) First solve Ly = b
(b) then solve U x = y.
Both systems are triangular and therefore can be easily solved using backward and forward
substitution.
We can carry out these two steps using the “” command in MATLAB:
y = Lb
x = Uy
However, in general, as observed before, MATLAB will return the LU factorization of P A
rather than A, thus, instead of solving the system Ax = b we will solve the equivalent system
P Ax = P b
(3)
To solve (3) using steps (a) and (b) defined above, we only need to replace b with P b:
y = L(P*b)
x = Uy
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
5
EXERCISE 4
Enter the matrix A and the vector b

−9
 3
A=
 −6
−1
in MATLAB:

6 6 −8
3 4 −9 
,
5 3 −6 
2 6
4


−85
 23 

b=
 −49 
−69
The exact solution to the system Ax = b is the vector x = (7, −2, −7, −4)T .
(a) Enter [L,U,P] = lu(A) to find the LU factorization of the matrix P A. Verify that
P A = LU by computing P A − LU .
(b) Use the LU factorization you found in part (a) to solve the system Ax = b. Call the
computed solution x lu (don’t forget that you will need P*b).
(c) Enter the vector x and compare your solution x lu from part (b) with the exact solution
x by computing norm(x lu – x)
NOTE: the norm function gives the magnitude of the vector,
p that is, for a vector a =
(a1 , a2 , . . . , an )T , the norm of a is defined as: norm(a) = a21 + a22 + . . . + a2n . Thus
the command norm(x lu – x) measures how close the computed solution x lu is to the
exact solution x. You should expect a very small number.
EXERCISE 5
In this exercise we will compare in MATLAB the speed of two methods for solving the system
Ax = b when A is an invertible square matrix: computing the RREF and computing the LU
factorization. Note that, although the number of operations for obtaining the LU factorization
is the same as the gaussian elimination, the LU factorization has the advantage that once the
matrix A is decomposed, the substitution step can be carried out efficiently for different values
of b. Thus the LU factorization is certainly preferable when solving the system Ax = b with
different values of b for the same A.
We will use the MATLAB tic and toc command to measure the computation times.
Enter:
A = rand(500);
x = ones(500,1);
b=A*x;
Important: Be sure to use semicolon after each command so that matrices and vectors are
not displayed. Do not print or include these large matrices and vectors in your lab write-up.
(a) Solve Ax = b using the reduced row echelon form and store the solution in x rref:
tic; R = rref([A, b]);
x_rref = R(:,end); toc
(make sure that all the commands are on the same line).
(b) Solve Ax = b using the LU factorization as you did in EXERCISE 4(a)(b), and calculate
the elapsed time using the tic toc commands. Store the solution in x lu.
Which method is faster?
NOTE: Make sure you use semicolon; the only output should be the elapsed time. Also,
don’t forget that, when using tic toc all the commands should be in one line, so you
should find L,U,P, solve for y, and solve for x lu all on one line.
(c) Compare the solutions from parts (a) and (b) with the exact solution x by computing
norm(x rref – x) and norm(x lu – x). How accurate are the solutions from parts (a)
and (b)?
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v1 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
6
MAT 275 Laboratory 3
Numerical Solutions by Euler and Improved Euler
Methods
In this session we look at basic numerical methods to help us understand the fundamentals of numerical
approximations. Our objective is as follows.
1. Implement Euler’s method as well as an improved version to numerically solve an IVP.
2. Compare the accuracy and efficiency of the methods with methods readily available in MATLAB.
3. Apply the methods to specific problems and investigate potential pitfalls of the methods.
Instructions: For your lab write-up, follow the instructions of LAB 1.
Euler’s Method
Consider the initial value problem (IVP): y 0 = f (t, y), y(t0 ) = y0 , where t0 , y0 are given numbers and f
is a given function of two variables. We are looking for the solution (i.e. function) y(t) on an interval
[t0 , tfinal ]. In many practical applications, the function f can be so complicated that it is not possible
to solve the IVP analytically. In such a situation we want to find a numerical solution. That is, for a
given number of points t0 > clear t y % no comma b / w t and y ! type ” help clear ” for more info
>> y (1) = -2; t (1) =0; h =0.25;
Here, the initialization y(1)=-2 should not be interpreted as “the value of y at 1 is −2”, but rather
as “the first value in the array y is −2”. In other words, the 1 in y(1) is an index, not a time value!
Unfortunately, MATLAB indices in arrays must be positive (a legacy from Fortran…).
Since f is simple enough we may use the anonymous function syntax:
>> f=@(t,y) 1.5*y;
The successive approximations at increasing values of t are then obtained as follows:
>> y (2)= y (1)+ h * f ( t (1) , y (1)) , t (2)= t (1)+ h % Euler approximation at t = 0.25
y =
-2.0000
-2.7500
t =
0
0.2500
>> y (3)= y (2)+ h * f ( t (2) , y (2)) , t (3)= t (2)+ h % Euler approximation at t = 0.5
y =
-2.0000
-2.7500
-3.7812
t =
0
0.2500
0.5000
The arrays y and t are now 1 × 3 row arrays. The dimension increases as new values of y and t are
computed.
It is obviously better to implement these steps in a for loop. The following loop implements the first 5
steps with the stepsize h = 0.25.
y (1) = -2; t (1) =0; h =0.25;
for n = 1:5
y ( n +1) = y ( n ) + h * f ( t ( n ) ,y ( n ) ) ;
t ( n +1) = t ( n ) + h ;
end
Note that the output in each command has been suppressed (with a ;). The list of computed y values vs
t values can be output at the end by typing:
>> [ t (:) ,y (:) ]
% same as [t ‘ ,y ‘]
ans =
0
-2.0000
0.2500
-2.7500
0.5000
-3.7812
0.7500
-5.1992
1.0000
-7.1489
1.2500
-9.8298
The next step is to write a function file (i.e. m-file) that takes in the following inputs: the name of the
function (f ) defining the ODE, the time span (i.e. the interval [t0 , tfinal ]), the initial condition (y0 ), and
THIS CONTENT IS PROTECTED AND MAY NOT BE SHARED, UPLOADED, SOLD OR DISTRIBUTED
2019 v6 Copyright@ School of Mathematical and Statistical Sciences, Arizona State University
2
the number of steps (N ) used in the approximation. Note that, if we are given the number of steps N ,
tfinal − t0
then the step size h can be easily computed as h =
.
N
As an output, we want the array of (t, y) points, where for each time step t, the corresponding y value is
the numerical approximation of the solution of the IVP: y 0 = f (t, y), y(t0 ) = y0 . The following Matlab
function (given in the file euler.m) implements these ideas.
function [t , y ] = euler (f , tspan , y0 , N )
%
%
%
%
%
%
%
%
%
%
%
%
%
Solves the IVP y ‘ = f (t , y ) , y ( t0 ) = y0 in the interval tspan =[ t0 , tf ]
using Euler ‘ s method with N time steps .
Input :
f = name of the anonymous function or M – file that evaluates
the ODE ( if not an anonymous function , use : euler (@ f , tspan , y0 , N ) )
For a system , the f must be given as column vector .
tspan = [ t0 , tf ] where t0 = initial time and tf = final time
y0 = initial value of the dependent variable . If solving a system ,
initial conditions must be given as a vector .
N = number of steps used .
Output :
t = vector of time values where the solution was computed
y = vector of computed solution values .
m = length ( y0 ) ;
t0 = tspan (1) ;
tf = tspan (2) ;
h = ( tf – t0 ) / N ;
% evaluate the time step size
t = linspace ( t0 , tf , N +1) ; % create the vector of t – values
y = zeros (m , N +1) ;
% allocate memory for the output y
y (: ,1) = y0 ‘;
% set initial condition
for n =1: N
y (: , n +1) = y (: , n ) + h * f ( t ( n ) ,y (: , n ) ) ; % implement Euler ‘ s method
end
t = t ‘; y = y ‘;
% change t and y from row to column vectors
end % end of the function euler
Remark: You should notice that the code above is slightly different from the first one we wrote (in
particular, note the use of “:” when creating the output y). Although the two codes are equivalent in
the scalar case, only the second one will work also for systems of Differential Equations (to be discussed
in Lab 4, when we will use this code again).
We can implement the function with, say, N = 70 steps in the interval [0, 1.75], by typing the following:
>> [t , y ] = euler (f ,[0 ,1.75] , -2 ,70) % use @ f if defined in separate m – file
>> [ t (1:5) , y (1:5)]
% print first five rows of the mat…
Purchase answer to see full
attachment