Question

4. (Applying LU and LUP decompositions) In this problem, we'll use the LU/LUP decomposition to solve a linear system of equations.

a) For A = [12 -8 13 -1 13;14 11 -5 -5 -7;1 -8 -9 10 8;-11 10 -8 3 8;-11 -8 4 2 -4] find matrices P, L, and U so that PA = LU using Matlab's lu function. Based on your results: did Matlab use pivoting during the lu-computation?

b) For b = [4;-4;-5;3;7] solve Ax = b using the LU decomposition as follows. Solving Ax = b is the same as solving PAx = Pb. (With P from a). Since PA = LU, we need to solve LUx = Pb, and we can split that into two triangular systems as follows: Ly = Pb, and Ux = y. Solve both of these systems using Matlab's linsolve, state x and y explicitly.

c) Compare the quality of the x you found in b to the solution of Ax = b you get from using linsolve. (As in 2d, work with the differences Ax - b).

d) You want so solve Ax = b for various vectors b, so you collect them into a single matrix B. So your goal is to find a matrix X so that AX = B (one column in X for each column in B). Working with P, L, and U from parts a/b we see that this amounts to solving two systems: LY = PB and UX = Y. For B = [18 -7 -14 10 -14 -2 13 12 -15 -15;-15 -14 4 -2 13 -16 15 -3 -15 14;3 12 -10 -17 2 19 -17 17 15 5;-1 -8 6 -11 20 -20 -4 -13 3 -6;-20 1 8 17 -17 11 -10 -10 2 1] solve these two equations using Matlab's linsolve.

First find Y in LY = PB, and then use that to find X in UX = Y. Check that AX - B is close to the zero matrix.

For this last problem, work with format short (or even format compact) so that the matrices don't use up too much screenspace.

Answer 1

%%Matlab code for solving linear system
clear all
close all

%A matrix
A = [12 -8 13 -1 13;14 11 -5 -5 -7;1 -8 -9 10 8;-11 10 -8 3 8;-11 -8 4 2 -4];
fprintf('A matrix is \n')
disp(vpa(A,3))

%LU of A
[L,U,P] = lu(A);
fprintf('L*U for A= ')
disp(vpa(L*U,3))
fprintf('Yes partial pivoting done for LU decomposition.\n')

%b vector
b = [4;-4;-5;3;7];
fprintf('b vector is \n')
disp(b)

%Ly=Pb
%Solution using linsolve
y=linsolve(L,P*b);
fprintf('Solution using linsolve for Ly=Pb for y= \n')
disp(y)
%Ux=y
%Solution using linsolve
x=linsolve(U,y);
fprintf('Solution using linsolve for Ux=y for x= \n')
disp(x)

%solution using linsolve for Ax=b
%Solution using linsolve
xx=linsolve(A,b);
fprintf('Solution using linsolve for Ax=b for xx= \n')
disp(xx)

%error in Ax=b and LU
err1=A*x-b;
fprintf('Error in LU solution\n')
disp(err1)

err2=A*xx-b;
fprintf('Error in Ax=b solution\n')
disp(err2)

%B matrix
B = [18 -7 -14 10 -14 -2 13 12 -15 -15;...
    -15 -14 4 -2 13 -16 15 -3 -15 14;...
    3 12 -10 -17 2 19 -17 17 15 5;...
    -1 -8 6 -11 20 -20 -4 -13 3 -6;...
    -20 1 8 17 -17 11 -10 -10 2 1] ;
fprintf('B matrix is\n')
disp(vpa(B,3))
%Ly=Pb
%Solution using linsolve
Y=linsolve(L,P*B);
fprintf('Solution using linsolve for LY=PB for Y= \n')
disp(vpa(Y,3))
%Ux=y
%Solution using linsolve
X=linsolve(U,Y);
fprintf('Solution using linsolve for UX=Y for X= \n')
disp(vpa(X,3))

%solution using linsolve for Ax=b
%Solution using linsolve
XX=linsolve(A,B);
fprintf('Solution using linsolve for AX=B for XX= \n')
disp(vpa(XX,3))

%error in Ax=b and LU
err1=A*X-B;
fprintf('Error in LU solution\n')
disp(vpa(err1,3))

%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%