Introduced forward_substitution()

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

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
+

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
+end

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)
+[U, Mk] = outer_product_gaussian_elimination(B);
+back_substitution(U, b);
+
+A = gaussian_elimination_with_complete_pivoting(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