Crout's LU Decomposition

This code is designed to solve a system of linear equations using Crouts LU decompostion. It seems to fall apart with some scenerios while "working" in others. Upper and Lower are correct, and I also believe Y is correct. Can someone verify if the issue lies within the U * Y = X portion of this code, or is it L * Y = B. Looking online I could not find a completed example.

A = [1 3 4 8

2 1 2 3

4 3 5 8

9 2 7 4];

B = [1

1

1

1];

matrixSize=length(A);

Lower=zeros(size(A));

Upper=zeros(size(A));

Lower(:,1)= A(:,1); %Set the first column of L to the frist column of A

Upper(1,:)=A(1,:)/Lower(1,1); % Create the first row of upper, divide by L(1,1)

Upper(1,1)=1; % Start the identity matrix

for k=2:matrixSize

for j=2:matrixSize

for i=j:matrixSize

Lower(i,j)=A(i,j) - dot(Lower(i,1:j-1),Upper(1:j-1,j));

end

Upper(k,j)=(A(k,j)-dot(Lower(k,1:k-1),Upper(1:k-1,j)))/Lower(k,k);

end

end

Upper

Lower

% L * Y = B

Y = zeros(matrixSize, 1);

% BASE CASE, SOLVE THE FIRST ONE

Y(1) = B(1);

for row = 2 : matrixSize %2 - number or rows

Y(row) = B(row);

for col = 1 : row - 1 %1 - row number

Y(row) = Y(row) - Lower(row, col) * Y(col);

end

Y(row) = Y(row) / Lower(row, row)

end

Y

% U * X = Y

X = zeros(matrixSize, 1);

X(matrixSize) = Y(matrixSize) / Upper(matrixSize,matrixSize);

for row = matrixSize - 1 : -1 : 1 %second to last row - 1

temp = 0;

for col = matrixSize : -1 : 1 %number of rows to row

temp = temp + Upper(row,col) * X(col);

end

X(row) = (Y(row) - temp) / Upper(row,row);

end

X

Answer 1

[Alar Kööbi]

% BASE CASE, SOLVE THE FIRST ONE

Y(1) = B(1)/Lower(1,1);%<------This one