I\'m trying to write a UDF in Matlab that will give me L and U for the LU decomp
ID: 3573273 • Letter: I
Question
I'm trying to write a UDF in Matlab that will give me L and U for the LU decomposition of the matrix A with the ones on the diagonal of the L matrix.
So far I have this code, but I keep getting NaN and INF answers in my matrixes for all of the sample ones that I'm trying.
Can you find where I'm going wrong on my thought process?
function [L, U] = metzler3_get_LU(A)
% The function decomposes the matrix A into a lower triangular matrix L
% and an upper triangular matrix U, using Crout's method such that A=LU.
% Input variable:
% A The matrix of coefficients.
% Output variables:
% L Lower triangular matrix.
% U Upper triangular matrix.
[R, C] = size(A);
L = zeros(R);
U = zeros(C);
if R ~= C
error('A is not a square matrix.');
end
if abs(det(A)) <= 10^-14;
error('Matrix is singular; the LU decomposition cannot be computed by this function.');
end
for j = 1:R
U(1, j) = A(1, j);
L(j, j) = 1;
end
for i = 2:R
L(i,1) = A(i, 1) / A(1, 1);
end
for j = 2:R-1
for i = 3:j
L(i, j) = (A(i, j) - L(i, 1:j - 1) * U(1:j - 1, j)) / U(j, j);
end
end
for j = 2:R
for i= j:R
U(i, j) = A(i, j) - L(i, 1:i - 1) * U(1:i - 1, j);
end
end
end % function [L, U] = metzler3_get_LU(A)
Explanation / Answer
Hii there,
Check this code .It may help you
matlab code
function[L R x]=LR2(A,b)
% This program will find a solution to Ax=b using first giving the decomposition of the matrix into L and R and then solving.
% Part 1 - Is this matrix square and nonsingular? give an error if not
[z y]=size(A);
if (z ~= y )
disp ( 'LR2 error: Matrix must be square' );
return;
end;
if det(A)==0
end
% Part 2 : Decomposition of matrix into L and R
L=zeros(z,z);
R=zeros(z,z);
for i=1:z
% Finding L
for k=1:i-1
L(i,k)=A(i,k);
for j=1:k-1
L(i,k)= L(i,k)-L(i,j)*R(j,k);
end
L(i,k) = L(i,k)/R(k,k);
end
% Finding R
for k=i:z
R(i,k) = A(i,k);
for j=1:i-1
R(i,k)= R(i,k)-L(i,j)*R(j,k);
end
end
end
for i=1:z
end
% Program shows R and L
R
L
% Now use a vector y to solve 'Ly=b'
y=zeros(z,1);
y(1)=b(1)/L(1,1);
for i=2:z
end;
% Now we use this y to solve Rx = y
x=zeros(z,1);
x(1)=y(1)/R(1,1);
for i=2:z
%print x x end