From 19cba61d44a69cee5edabf2329a74381c30eae41 Mon Sep 17 00:00:00 2001 From: EdwardEisenhauer Date: Sat, 13 Mar 2021 20:23:23 +0100 Subject: [PATCH] Added Algorithm 8 and refactored code --- .gitignore | 1 + .../Alg1_outer_product_gaussian_elimination.m | 17 +++--- ...ssian_elimination_with_complete_pivoting.m | 18 +++---- .../Alg3_forward_substitution.m | 27 +++++----- .../Alg4_back_substitution.m | 2 +- .../Alg5_gauss_jordan_elimination.m | 18 +++---- .../Alg6_RREF.m | 4 +- .../Alg7.m | 25 ++++----- .../Alg8.m | 54 +++++++++++++++++++ 9 files changed, 110 insertions(+), 56 deletions(-) create mode 100644 Direct Methods for Solving Linear Systems/Alg8.m diff --git a/.gitignore b/.gitignore index 3e3461a..961f132 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ +# MATLAB *.asv diff --git a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m index 565a867..31648c1 100644 --- a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m +++ b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m @@ -1,20 +1,21 @@ function A = Alg1_outer_product_gaussian_elimination(A) -% ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1) -% Matrix A has to be a square singular matrix. +% Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1) +% Performs a gaussian eliminaion on a square matrix A. -[m, n] = size(A); -if m ~= n - error('Matrix is not squared!') -end + [m, n] = size(A); + + if m ~= n + error('Matrix is not squared!') + end % if det(A) == 0 % error('Matrix is not nonsingular!') % end - for k = 1 : m-1 + for k = 1 : m-1 rows = k + 1 : m; A(rows, k) = A(rows, k)/A(k, k); A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows); - end + end end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m index f48348b..7c01c8e 100644 --- a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m +++ b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m @@ -1,27 +1,27 @@ function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) -% ADDME Algorithm 2: Gaussian Elimination with Complete Pivoting. +% Algorithm 2: Gaussian Elimination with Complete Pivoting. % [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) [n, m] = size(A); + if n ~= m error('Matrix is not squared!') end -if det(A) == 0 - error('Matrix is not nonsingular!') -end +% if det(A) == 0 +% error('Matrix is not nonsingular!') +% end p = 1:n; q = 1:n; -% for k = 1 : n-1 for k = 1 : n-1 i = k:n; j = k:n; [max_val, rows_of_max_in_col] = max(abs(A(i, j))); - [max_val, max_col] = max(max_val); + [~, max_col] = max(max_val); max_row = rows_of_max_in_col(max_col); - % Assing value of mi and lambda in respect to the main A matrix + % Assign value of mu and lambda in respect to the main A matrix [mi, lm] = deal(max_row+k-1, max_col+k-1); A([k mi], 1:n) = deal(A([mi k], 1:n)); A(1:n, [k lm]) = deal(A(1:n, [lm k])); @@ -36,8 +36,8 @@ for k = 1 : n-1 end end -U = triu(A); -L = tril(A, -1) + eye(n); I = eye(n); +U = triu(A); +L = tril(A, -1) + I; P = I(p, :); Q = I(:, q); \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m index 84a0acb..472da08 100644 --- a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m +++ b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m @@ -1,17 +1,18 @@ -% Algorithm 3: Forward Substitution (Alg. 3.1.1) function b = Alg3_forward_substitution(L, b) +% Algorithm 3: Forward Substitution (Golub, Loan, Alg. 3.1.1) -[n, m] = size(L); -if n ~= m - error('Matrix is not squared!') -end - -if length(b) ~= n - error('Vector b has wrong length!') -end - -b(1, :) = b(1, :)/L(1,1); -for i = 2:n - b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i); + [m, n] = size(L); + if m ~= n + error('Matrix is not squared!') + end + + if length(b) ~= m + error('Vector b has wrong length!') + end + + b(1, :) = b(1, :)/L(1,1); + for i = 2:m + b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i); + end end diff --git a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m index 622afe1..9f25ebf 100644 --- a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m +++ b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m @@ -1,5 +1,5 @@ function b = Alg4_back_substitution(U,b) -% ADDME Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2) +% Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2) % Returns vetor b with solution to he Ux = b. [m, n] = size(U); diff --git a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m index a48b804..f510427 100644 --- a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m +++ b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m @@ -2,19 +2,17 @@ function A = Alg5_gauss_jordan_elimination(A) % Algorithm 5: Gauss-Jordan Elimination % Argument A is an augmented matrix -% M – rows, N – columns +[m, n] = size(A); -[M, N] = size(A); - -for m = 1 : M +for k = 1 : m - row = A(m, :); - row = row/row(m); - A(m, :) = row; + row = A(k, :); + row = row/row(k); + A(k, :) = row; - for n = 1 : M - if n ~= m - A(n, :) = A(n, :)-(A(n, m))*row; + for i = 1 : m + if i ~= k + A(i, :) = A(i, :)-(A(i, k))*row; end end end diff --git a/Direct Methods for Solving Linear Systems/Alg6_RREF.m b/Direct Methods for Solving Linear Systems/Alg6_RREF.m index d2cec19..5d055db 100644 --- a/Direct Methods for Solving Linear Systems/Alg6_RREF.m +++ b/Direct Methods for Solving Linear Systems/Alg6_RREF.m @@ -1,9 +1,7 @@ function A = Alg6_RREF(A) -% ADDME Algorithm 6: Reduced Row Echelon Form (RREF) +% Algorithm 6: Reduced Row Echelon Form (RREF) % A = Alg6_RREF(A) returns RREF of matrix A. -% M – rows, N – columns - [M, N] = size(A); n = 0; diff --git a/Direct Methods for Solving Linear Systems/Alg7.m b/Direct Methods for Solving Linear Systems/Alg7.m index 6e8e9ca..40f4726 100644 --- a/Direct Methods for Solving Linear Systems/Alg7.m +++ b/Direct Methods for Solving Linear Systems/Alg7.m @@ -1,7 +1,8 @@ function [L, U] = Alg7(A) -% ADDME LU Factorization without pivoting -% [L, U] = Alg7(A) decomposes matrix A into U – upper triangular matrix and -% L – lower unit triangular matrix such, that A = LU. +% LU Factorization without pivoting +% [L, U] = Alg7(A) decomposes matrix A using Crout's algorithm into +% U – upper triangular matrix and L – unit lower triangular matrix +% such that A = LU. [m, n] = size(A); @@ -14,21 +15,21 @@ function [L, U] = Alg7(A) % L = tril(A, -1) + eye(m); % Instead, we should use the Crout's algorithm -% It returns an Unit Upper Triangular matrix and a Lower Traingular matrix -% (in oppose to the Gaussian Elimination). -A -L = zeros(m, n) -U = eye(m, n) +L = zeros(m, n); +U = eye(m, n); for k = 1:n - for i = k:n - L(i, k) = A(i, k) - dot(L(i,1:k-1), U(1:k-1, k)); - end for j = k : n - U(k, j) = (A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j))) / L(k, k); + U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j)); end + + for i = k:n + L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k); + end + + end end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg8.m b/Direct Methods for Solving Linear Systems/Alg8.m new file mode 100644 index 0000000..82666be --- /dev/null +++ b/Direct Methods for Solving Linear Systems/Alg8.m @@ -0,0 +1,54 @@ +function [L, U, P] = Alg8(A) +% LU Factorization with partial pivoting +% [L, U, P] = Alg8(A) decomposes matrix A using Crout's algorithm into +% U – upper triangular matrix and L – unit lower triangular matrix +% using partial pivoting, interchanging A's rows, such that AP = LU. + +[m, n] = size(A); + +L = zeros(m, n); +U = eye(m, n); +P = eye(m, n); + +for k = 1:m + if A(k, k) == 0 + non_zero_col = find(A(k, :), 1); + A(:, [k non_zero_col]) = deal(A(:, [non_zero_col k])); + P(:, [k non_zero_col]) = deal(P(:, [non_zero_col k])); + end +end + +% TODO: Permutation A*P does not work + + +% +% for i = m + 1 : M +% A(i:end, m+1:end); % Partial matrix (in which we are looking for non-zero pivots) +% A(i:end, m+1); % Left-most column +% if ~any(A(i:end, m+1)) % If the left-most column has only zeros check the next one +% m = m + 1; +% end +% A(i:end, m+1:end); +% if A(i, m+1) == 0 +% non_zero_row = find(A(i:end,m+1), 1); +% if isempty(non_zero_row) +% continue +% end +% A([i, i+non_zero_row-1], :) = deal(A([i+non_zero_row-1, i], :)); +% end +% end + +for k = 1:n + + for j = k : n + U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j)); + end + + for i = k:n + L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k); + end + + +end + +end \ No newline at end of file