Develop an Mfile to determine LU factorization of a square m

Develop an M-file to determine LU factorization of a square matrix without partial pivoting. That is, develop a function that is passed the square matrix and returns the triangular matrices [L] and [U]. Test your function by using it to solve a system of equations as follows: 10x_1 + 2x_2 - x_3 = 27 -3x_1 - 6x_2 + 2x_3 = -61.5 x_1 + x_2 + 5x_3 = -21.5 Confirm that your function is working properly by verifying that [L] [U] = [A] and by using the built-in function lu.

Solution

% LU % A*x=b -> A = L*U % L*U*x = b -> U*x = y % L*y = b clear,clc %Input definitions A = [ 10 2 -1; -3 -6 2; 1 1 5]; % Square matrix \"A\" -> A*x=b b = [27; -61.5; -21.5]; % Independente vector \"b\" n = length(A); U = [A]; N = 1:n; L = diag(ones(1,n)); %Pivoting... commentted because is not wanted % for k = 1:n-1 % for p = k+1:n % pivote = U(k,k); % if abs(U(p,k)) > abs(pivote) % for j = 1:n % next = U(k,j); % U(k,j) = U(p,j); % U(p,j) = next; % end % nextb = b(k); % b(k) = b(p); % b(p) = nextb; % end % end % end %Factorization LU for k = 1:n-1 for i = k+1:n factor = U(i,k)/U(k,k); L(i,k) = factor; for j = k:n U(i,j) = U(i,j) - factor*U(k,j); end end end %Front-Sustitution y(1) = b(1); for i = 2:n sum = b(i); for j = 1:i-1 sum = sum - L(i,j)*y(j); end y(i) = sum; end %Back-Sustitution x(n) = y(n)/U(n,n); for i = n-1:-1:1 sum = y(i); for j = i+1:n sum = sum - U(i,j)*x(j); end x(i) = (sum)/U(i,i); end %Resultados table = [N;x]; disp(\'The solution is\') fprintf(\'x%d = %9.4f \ \',table) L U disp(\'L*U\') L*U