From 87f957d90306746173854a61804f2a5c164eaa5a Mon Sep 17 00:00:00 2001 From: EdwardEisenhauer Date: Sat, 6 Mar 2021 18:17:36 +0100 Subject: [PATCH] Introduced forward_substitution() --- .../back_substitution.m | 6 ++--- .../forward_substitution.m | 17 ++++++++++++ ...ssian_elimination_with_complete_pivoting.m | 27 ++++++++++--------- .../main.m | 10 ++++--- .../outer_product_gaussian_elimination.m | 7 ++--- 5 files changed, 42 insertions(+), 25 deletions(-) create mode 100644 Direct Methods for Solving Linear Systems/forward_substitution.m diff --git a/Direct Methods for Solving Linear Systems/back_substitution.m b/Direct Methods for Solving Linear Systems/back_substitution.m index bbe1a07..42de1bf 100644 --- a/Direct Methods for Solving Linear Systems/back_substitution.m +++ b/Direct Methods for Solving Linear Systems/back_substitution.m @@ -1,5 +1,5 @@ % Argorithm 4: Back Substitution (Alg. 3.1.2) -function back_substitution(U,b) +function b = back_substitution(U,b) [n, m] = size(U); if n ~= m @@ -14,9 +14,9 @@ if det(U) == 0 error('Matrix is not nonsingular!') end -b(n) = b(n)/U(n, n) +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 new file mode 100644 index 0000000..2a01768 --- /dev/null +++ b/Direct Methods for Solving Linear Systems/forward_substitution.m @@ -0,0 +1,17 @@ +% Algorithm 3: Forward Substitution (Alg. 3.1.1) +function b = forward_substitution(L, b) + +[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); +end + 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 b835b52..91be0a6 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 gaussian_elimination_with_complete_pivoting(A) +function A = gaussian_elimination_with_complete_pivoting(A) % det(A(mi, lm)) @@ -8,23 +8,24 @@ if n ~= m end if det(A) == 0 - error('Matrix is not nonsingular!') + error('Matrix is not nonsingular!') end 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, 1:n) = A(mi, 1:n) - A(1:n, k) = A(1:n, lm) - p(k) = mi - q(k) = lm + 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) + 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 -end \ No newline at end of file +end diff --git a/Direct Methods for Solving Linear Systems/main.m b/Direct Methods for Solving Linear Systems/main.m index 37da31a..7104b5c 100644 --- a/Direct Methods for Solving Linear Systems/main.m +++ b/Direct Methods for Solving Linear Systems/main.m @@ -1,8 +1,10 @@ clear all; -B = [2, -1, 0, 0; -1, 2, -1, 0;0, -1, 2, -1; 0, 0, -1, 2] +B = [2, -1, 0, 0; -1, 2, -1, 0;0, -1, 2, -1; 0, 0, -1, 2]; -b = [0;0;0;5] +b = [0;0;0;5]; -[U, Mk] = outer_product_gaussian_elimination(B) -back_substitution(U, b) \ No newline at end of file +[U, Mk] = outer_product_gaussian_elimination(B); +back_substitution(U, b); + +A = gaussian_elimination_with_complete_pivoting(B); 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 83a811b..085a4c8 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 @@ -10,14 +10,11 @@ end error('Matrix is not nonsingular!') end - A; - for k = 1 : n-1 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); - A; end -U = triu(A) -Mk = diag(A,-1) % Gauss vector +U = triu(A); +Mk = diag(A,-1); % Gauss vector