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 64dd6f7..1e47f17 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,14 +1,15 @@ -% Algorithm 1: Outer Product Gaussian Elimination (Alg. 3.2.1) function A = Alg1_outer_product_gaussian_elimination(A) +% ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1) + [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 for k = 1 : n-1 rows = k + 1 : n; 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 172bc67..f48348b 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,4 +1,6 @@ function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) +% ADDME 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 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 f994f62..622afe1 100644 --- a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m +++ b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m @@ -1,24 +1,30 @@ -% Argorithm 4: Back Substitution (Alg. 3.1.2) function b = Alg4_back_substitution(U,b) +% ADDME Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2) +% Returns vetor b with solution to he Ux = b. -[n, m] = size(U); -if n ~= m +[m, n] = size(U); + +if U ~= triu(U) + error('Matrix is not upper triangular!') +end + +if m ~= n error('Matrix is not squared!') end -if length(b) ~= n +if length(b) ~= m error('Vector b has wrong length!') end -if det(U) == 0 - error('Matrix is not nonsingular!') -end +% if det(U) < 0.001 +% error('Matrix is not nonsingular!') +% end -% b(n, :) so that matrices are also accepted +% b(m, :) so that matrices are also accepted -b(n, :) = b(n, :)/U(n, n); -for i = n-1:-1:1 - b(i, :) = (b(i, :) - U(i, i+1 : n)*b(i+1 : n, :))/U(i, i); +b(m, :) = b(m, :)/U(m, m); +for i = m-1:-1:1 + b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i); end end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg6_RREF.m b/Direct Methods for Solving Linear Systems/Alg6_RREF.m index 01314d1..d2cec19 100644 --- a/Direct Methods for Solving Linear Systems/Alg6_RREF.m +++ b/Direct Methods for Solving Linear Systems/Alg6_RREF.m @@ -1,5 +1,6 @@ function A = Alg6_RREF(A) -% Algorithm 6: Reduced Row Echelon Form (RREF) +% ADDME Algorithm 6: Reduced Row Echelon Form (RREF) +% A = Alg6_RREF(A) returns RREF of matrix A. % M – rows, N – columns @@ -8,17 +9,15 @@ function A = Alg6_RREF(A) n = 0; for m = 1 : M - n = n + 1 + n = n + 1; if n > N break end - A % We want the left-most coefficient to be 1 (pivot) row = A(m, :); if row(m) == 0 - n = n + 1 + n = n + 1; end - [m ,n] row = row/row(n); A(m, :) = row; @@ -28,8 +27,6 @@ for m = 1 : M end end - A - 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 diff --git a/Direct Methods for Solving Linear Systems/Alg7.m b/Direct Methods for Solving Linear Systems/Alg7.m new file mode 100644 index 0000000..ecb620a --- /dev/null +++ b/Direct Methods for Solving Linear Systems/Alg7.m @@ -0,0 +1,13 @@ +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. + +[m, n] = size(A); + +A = Alg1_outer_product_gaussian_elimination(A); + +U = triu(A); +L = tril(A, -1) + eye(m); + +end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/main.m b/Direct Methods for Solving Linear Systems/main.m index 243cd7f..2ae813a 100644 --- a/Direct Methods for Solving Linear Systems/main.m +++ b/Direct Methods for Solving Linear Systems/main.m @@ -1,4 +1,5 @@ clear all; +clc; A = [2, -1, 0, 0; -1, 2, -1, 0; @@ -9,8 +10,7 @@ b = [0;0;0;5]; B = Alg1_outer_product_gaussian_elimination(A); U = triu(B); -L = tril(B, -1); -x = back_substitution(U, b); +x = Alg4_back_substitution(U, b); %% Problem 1 clear all; @@ -44,9 +44,9 @@ b = [1;2;1]; b = P*b; % Ly = b and Ux = y -y = forward_substitution(L, b); +y = Alg3_forward_substitution(L, b); -x = Q*back_substitution(U, y) +x = Q*Alg4_back_substitution(U, y) % L*U