From 82e6cf7d1046d6cee16f7e8b044ec33e7ec6c4b7 Mon Sep 17 00:00:00 2001
From: Nguyễn Gia Phong <vn.mcsinyx@gmail.com>
Date: Sun, 16 Feb 2020 14:26:55 +0700
Subject: [usth] Numerical Method is MATH2.4

---
 usth/MATH2.4/labwork/3/GaussPivot.m | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 usth/MATH2.4/labwork/3/GaussPivot.m

(limited to 'usth/MATH2.4/labwork/3/GaussPivot.m')

diff --git a/usth/MATH2.4/labwork/3/GaussPivot.m b/usth/MATH2.4/labwork/3/GaussPivot.m
new file mode 100644
index 0000000..7764b3a
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/GaussPivot.m
@@ -0,0 +1,25 @@
+function x = GaussPivot (A, b)
+  [m n] = size (A);
+  m ~= n && error ('Matrix A must be square');
+  nb = n + 1;
+  Aug = [A b];
+  % forward elimination
+  for k = 1 : n - 1
+    % partial pivoting
+    [big i] = max (abs (Aug(k:n, k)));
+    ipr = i + k - 1;
+    if ipr ~= k
+      Aug([k ipr], :) = Aug([ipr k], :);
+    end
+    for i = k + 1 : n
+      coeff = Aug(i, k) / Aug(k, k);
+      Aug(i, k:nb) = Aug(i, k:nb) - coeff * Aug(k, k:nb);
+    end
+  end
+  % back substitution
+  x = zeros (n, 1);
+  x(n) = Aug(n, nb) / Aug(n, n);
+  for i = n - 1 : -1 : 1
+    x(i) = (Aug(i, nb) - Aug(i, i+1:n) * x(i+1:n)) / Aug(i, i);
+  end
+end
-- 
cgit 1.4.1