From d4d323fd19b28cc0b01de055fd11ff6c783b55d7 Mon Sep 17 00:00:00 2001 From: Sergiusz Warga Date: Sun, 7 Mar 2021 22:56:18 +0100 Subject: [PATCH] All algorithms from 1 to 5 work just fine! --- .../back_substitution.m | 6 +- .../forward_substitution.m | 4 +- .../gauss_jordan_elimination.m | 38 ++++----- ...ssian_elimination_with_complete_pivoting.m | 25 +++--- .../main.m | 80 +++++++++++++++++-- .../outer_product_gaussian_elimination.m | 7 +- 6 files changed, 108 insertions(+), 52 deletions(-) diff --git a/Direct Methods for Solving Linear Systems/back_substitution.m b/Direct Methods for Solving Linear Systems/back_substitution.m index 42de1bf..0861344 100644 --- a/Direct Methods for Solving Linear Systems/back_substitution.m +++ b/Direct Methods for Solving Linear Systems/back_substitution.m @@ -14,9 +14,11 @@ if det(U) == 0 error('Matrix is not nonsingular!') end -b(n) = b(n)/U(n, n); +% b(n, :) 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(i, :) = (b(i, :) - U(i, i+1 : n)*b(i+1 : n, :))/U(i, i); end end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/forward_substitution.m b/Direct Methods for Solving Linear Systems/forward_substitution.m index 2a01768..591e6d0 100644 --- a/Direct Methods for Solving Linear Systems/forward_substitution.m +++ b/Direct Methods for Solving Linear Systems/forward_substitution.m @@ -10,8 +10,8 @@ if length(b) ~= n error('Vector b has wrong length!') end -b(1) = b(1)/L(1,1); +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); + b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i); end diff --git a/Direct Methods for Solving Linear Systems/gauss_jordan_elimination.m b/Direct Methods for Solving Linear Systems/gauss_jordan_elimination.m index 20e8543..d7bf51c 100644 --- a/Direct Methods for Solving Linear Systems/gauss_jordan_elimination.m +++ b/Direct Methods for Solving Linear Systems/gauss_jordan_elimination.m @@ -1,31 +1,25 @@ +% Algorithm 5: Gauss-Jordan Elimination +% Input A is an augmented matrix function A = gauss_jordan_elimination(A) [n, m] = size(A); -if n ~= m + +if n + 1 ~= 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 -A - -for k = 1 : n-1 - i = k:n; - j = k:n; - A(i, j); - maximum = max(abs(A(i, j)), [], 'all'); - max_idx = find(abs(A==maximum)); - [mi, lm] = ind2sub(size(A), max_idx(1)); - A([k mi], 1:n) = deal(A([mi k], 1:n)); - A(1:n, [k lm]) = deal(A(1:n, [lm k])); - p(k) = mi; - q(k) = lm; - % Perform Gaussian elimination with the greatest pivot - if A(k, k) ~= 0 - rows = k+1 : n; - A(rows, k) = A(rows, k)/A(k, k); - A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows); +for k = 1 : m-1 + + row = A(k, :); + row = row/row(k); + A(k, :) = row; + for l = 1 : m-1 + if l ~= k + A(l, :) = A(l, :)-(A(l, k))*row; + end end end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/gaussian_elimination_with_complete_pivoting.m b/Direct Methods for Solving Linear Systems/gaussian_elimination_with_complete_pivoting.m index d566496..521d913 100644 --- a/Direct Methods for Solving Linear Systems/gaussian_elimination_with_complete_pivoting.m +++ b/Direct Methods for Solving Linear Systems/gaussian_elimination_with_complete_pivoting.m @@ -1,4 +1,4 @@ -function [U, L] = gaussian_elimination_with_complete_pivoting(A) +function [P, Q, L, U] = gaussian_elimination_with_complete_pivoting(A) [n, m] = size(A); if n ~= m @@ -9,24 +9,25 @@ if det(A) == 0 error('Matrix is not nonsingular!') end -A +p = 1:n; +q = 1:n; % for k = 1 : n-1 for k = 1 : n-1 i = k:n; j = k:n; - A(i, j) + A(i, j); [max_val, rows_of_max_in_col] = max(abs(A(i, j))); [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 - [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])) - p(k) = mi - q(k) = lm - % Perform Gaussian elimination with the greatest pivot + [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])); + p([k, mi]) = p([mi, k]); + q([k, lm]) = q([lm, k]); + % Perform Gaussian elimination with the greatest pivot if A(k, k) ~= 0 rows = k+1 : n; A(rows, k) = A(rows, k)/A(k, k); @@ -36,8 +37,6 @@ end U = triu(A); L = tril(A, -1) + eye(n); -p I = eye(n); -P = I(p, :) -q -Q = I(:, q) \ No newline at end of file +P = I(p, :); +Q = I(:, q); \ 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 9a23c13..129d201 100644 --- a/Direct Methods for Solving Linear Systems/main.m +++ b/Direct Methods for Solving Linear Systems/main.m @@ -1,15 +1,79 @@ clear all; -B = [2, -1, 0, 0; +A = [2, -1, 0, 0; -1, 2, -1, 0; - 3, -1, 2, -1; - 0, 4, -1, 2]; + 0, -1, 2, -1; + 0, 0, -1, 2]; b = [0;0;0;5]; -[U, L] = outer_product_gaussian_elimination(B); -back_substitution(U, b); +B = outer_product_gaussian_elimination(A); +U = triu(B); +L = tril(B, -1); +x = back_substitution(U, b); -[U, L] = gaussian_elimination_with_complete_pivoting(B) -L*U -% A = gauss_jordan_elimination(B) \ No newline at end of file +%% Problem 1 +clear all; +A = [2, -1, 0, 0; + -1, 2, -1, 0; + 0, -1, 2, -1; + 0, 0, -1, 2]; + +b = [0;0;0;5]; + +B = gauss_jordan_elimination([A b]) + +[P, Q, L, U] = gaussian_elimination_with_complete_pivoting(A); + +b = P*b; +% Ly = b and Ux = y +y = forward_substitution(L, b); + +x = Q*back_substitution(U, y); + +% L*U + +%% Problem 2 +A = [1, 1, 1; + 1, 1, 2; + 1, 2, 2]; + +b = [1;2;1]; + +[P, Q, L, U] = gaussian_elimination_with_complete_pivoting(A); + +b = P*b; +% Ly = b and Ux = y +y = forward_substitution(L, b); + +x = Q*back_substitution(U, y) + +% L*U + + +%% Problem 4 + +A = [0.835, 0.667; + 0.333, 0.266]; +b = [0.168; 0.067]; +bp = [0.168; 0.066]; + +kappa = cond(A) + +B = gauss_jordan_elimination([A b]) +Bp = gauss_jordan_elimination([A bp]) + +%% Problem 5 +% AX = I3 +A = [2, 1, 2; + 1, 2, 3; + 4, 1, 2]; + +[P, Q, L, U] = gaussian_elimination_with_complete_pivoting(A); + +I = P*eye(3); +% Ly = b and Ux = y +y = forward_substitution(L, I); + +X = Q*back_substitution(U, y) +inv(A) \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m b/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m index 6799cdd..8d43864 100644 --- a/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m +++ b/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m @@ -1,5 +1,5 @@ % Algorithm 1: Outer Product Gaussian Elimination (Alg. 3.2.1) -function [U, L] = outer_product_gaussian_elimination(A) +function A = outer_product_gaussian_elimination(A) [n, m] = size(A); if n ~= m @@ -14,7 +14,4 @@ end rows = k + 1 : n; A(rows, k) = A(rows, k)/A(k, k); A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows); - end - -U = triu(A); -L = tril(A, -1) + eye(n); + end \ No newline at end of file