cant find problem in my LU decomposition code

i wrote this code to determine my l and u for the given matrix,but when i calculate l*u some arrays in its last column isn't the same as given matrix,it has messed up my mind...i appreciate your helps to make my code work.

 if true
    % clear
clc
a=[0 3 -4 2;2 6 4 -3;-1 -1 2 3;3 0 0 -5];
c=[26;9;7;-19];
n=length(a);
for k=1:n;
    pivot=k;
    big=abs(a(k,k));
    for i=(k+1):n;
        if a(i,k)>big;
            big=a(i,k);
            pivot=i;
        end
    end
if pivot~=k;
    a([pivot,k],:)=a([k,pivot],:);
end
end
for i=1:n;
    l(i,1)=a(i,1);
end
for i=1:n;
    u(i,i)=1;
end
for j=2:n;
    u(1,j)=a(1,j)/l(1,1);
end
for j=2:(n-1);
    for i=j:n;
        sum=a(i,j);
        for k=1:(j-1);
        sum=sum-l(i,k)*u(k,j);
        end
        l(i,j)=sum;
    end

    for k=(j+1):n;
        sum1=a(j,k);
        for i=1:(j-1);
        sum1=sum1-l(j,i)*a(i,k);
        end
        u(j,k)=sum1/l(j,j);
    end
end
for j=n;
    sum2=a(j,j);
    for k=1:(n-1);
        sum2=sum2-l(j,k)*a(k,j);
    end
    l(j,j)=sum2;
end

Answer 1

[Naga]

Hello Hosein,

There are a couple of issues in your code that might cause incorrect results. You are checking 'a(i,k) > big', but you should be checking 'abs(a(i,k)) > big' for partial pivoting. Additionally, in the logic for computing (L) and (U), the computation of (U) should use the updated values of (L).

Here's a revised version of your code:

a = [0 3 -4 2; 2 6 4 -3; -1 -1 2 3; 3 0 0 -5];
n = length(a);

% Pivoting
for k = 1:n
    pivot = k;
    big = abs(a(k,k));
    for i = (k+1):n
        if abs(a(i,k)) > big
            big = abs(a(i,k));
            pivot = i;
        end
    end
    if pivot ~= k
        a([pivot,k],:) = a([k,pivot],:);
    end
end

% LU Decomposition
for i = 1:n
    L(i,i) = 1; % Diagonal of L is 1
    for j = i:n
        sum = 0;
        for k = 1:(i-1)
            sum = sum + L(i,k) * U(k,j);
        end
        U(i,j) = a(i,j) - sum;
    end
    for j = (i+1):n
        sum = 0;
        for k = 1:(i-1)
            sum = sum + L(j,k) * U(k,i);
        end
        L(j,i) = (a(j,i) - sum) / U(i,i);
    end
end