%% LU Decomposition % % This script constructs a random matrix A and random right hand side b, % finds the LU factorization with partial pivoting, and then solves using % forward and back-substitution. % function x = solve(A,b) [L,U] = my_lu(A); y = forward_solve(L,b); x = backward_solve(U,y); %% function [L,U] = my_lu(A) function [L,U] = my_lu(A) N = length(A); L = eye(N); U = A; for j = 1:N, for k = j+1:N, mp = U(k,j)/U(j,j); L(k,j) = mp; U(k,1:j) = 0; U(k,j+1:N) = U(k,j+1:N) - mp*U(j,j+1:N); end end end %% function [x] = forward_solve(L,U,b); function y = forward_solve(L,b) N = length(L); y = zeros(1,N); for i = 1:N, s = sum(L(i,1:i-1).*y(1:i-1)); y(i) = b(i) - s; end y = y(:); end function x = backward_solve(U,y) % Backward substitution - solve Ux = y N = length(y); x = zeros(1,N); for i = N:-1:1, s = sum(U(i,i+1:N).*x(i+1:N)); x(i) = (y(i) - s)/U(i,i); end x = x(:); end end