All algorithms from 1 to 5 work just fine!

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

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
 

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
+% if det(A) == 0
-    error('Matrix is not nonsingular!') 
+%     error('Matrix is not nonsingular!') 
-end
+% end
 
-A
+for k = 1 : m-1
     
-for k = 1 : n-1
+    row = A(k, :);
-    i = k:n;
+    row = row/row(k);
-    j = k:n;
+    A(k, :) = row;
-    A(i, j);
+    for l = 1 : m-1
-    maximum = max(abs(A(i, j)), [], 'all');
+        if l ~= k
-    max_idx = find(abs(A==maximum));
+            A(l, :) = A(l, :)-(A(l, k))*row;
-    [mi, lm] = ind2sub(size(A), max_idx(1));
+        end
-    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);
     end
 end

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)
+    [mi, lm] = deal(max_row+k-1, max_col+k-1);
-    A([k mi], 1:n) = deal(A([mi k], 1:n))
+    A([k mi], 1:n) = deal(A([mi k], 1:n));
-    A(1:n, [k lm]) = deal(A(1:n, [lm k]))
+    A(1:n, [k lm]) = deal(A(1:n, [lm k]));
-    p(k) = mi
+    p([k, mi]) = p([mi, k]);
-    q(k) = lm
+    q([k, lm]) = q([lm, k]);
-    % Perform Gaussian elimination with the greatest pivot
     
+    % 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, :)
+P = I(p, :);
-q
+Q = I(:, q);
-Q = I(:, q)

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, -1, 2, -1;
-    0, 4, -1, 2];
+    0, 0, -1, 2];
 
 b = [0;0;0;5];
 
-[U, L] = outer_product_gaussian_elimination(B);
+B = outer_product_gaussian_elimination(A);
-back_substitution(U, b);
+U = triu(B);
+L = tril(B, -1);
+x = back_substitution(U, b);
 
-[U, L] = gaussian_elimination_with_complete_pivoting(B)
+%% Problem 1
-L*U
+clear all;
-% A = gauss_jordan_elimination(B)
+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)

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
@@ -15,6 +15,3 @@ end
      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);