[usth] Numerical Method is MATH2.4

6 files changed, 84 insertions, 0 deletions

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
diff --git a/usth/MATH2.4/labwork/3/LU.m b/usth/MATH2.4/labwork/3/LU.m
new file mode 100644
index 0000000..860486d
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/LU.m
@@ -0,0 +1,10 @@
+function [L U] = LU (A)
+  [n n1] = size (A);
+  [L U] = deal (eye (n), A);
+  for k = 1:n
+    for i = k + 1 : n
+      L(i, k) = U(i, k) / U(k, k);
+      U(i, :) = U(i, :) - L(i, k)*U(k, :);
+    end
+  end
+end
diff --git a/usth/MATH2.4/labwork/3/LU_pivot.m b/usth/MATH2.4/labwork/3/LU_pivot.m
new file mode 100644
index 0000000..c506697
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/LU_pivot.m
@@ -0,0 +1,20 @@
+function [L U P] = LU_pivot (A)
+  [n _] = size (A);
+  [L P U] = deal (eye (n), eye (n), A);
+  for k = 1:n
+    [pivot m] = max (abs (U(k:n, k)));
+    m = m + k - 1;
+    if (m ~= k)
+      U([m k], :) = U([k m], :);  % interchange rows m and k in U
+      P([m k], :) = P([k m], :);  % interchange rows m and k in P
+      if k >= 2;  % very important point
+        % interchange rows m and k in columns 1:k-1 of L
+        L([m k], 1:k-1) = L([k m], 1:k-1);
+      end
+    end
+    for i = k + 1 : n
+      L(i, k) = U(i, k) / U(k, k);
+      U(i, :) -= L(i, k)*U(k, :);
+    end
+  end
+end
diff --git a/usth/MATH2.4/labwork/3/gauss.m b/usth/MATH2.4/labwork/3/gauss.m
new file mode 100644
index 0000000..09128dc
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/gauss.m
@@ -0,0 +1,12 @@
+function x = gauss (A, b)
+  issquare (A) || error ('Matrix A must be square');
+  Aug = [A b];
+  [m n] = size (Aug);
+  % forward elimination
+  for k = 1 : m - 1
+    for l = k + 1 : m
+      Aug(l, k:n) = Aug(l, k:n) - Aug(k, k:n)*Aug(l, k)/Aug(k, k);
+    end
+  end
+  x = ubst (Aug);
+end
diff --git a/usth/MATH2.4/labwork/3/labwork3.pdf b/usth/MATH2.4/labwork/3/labwork3.pdf
new file mode 100644
index 0000000..4e0c7e0
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/labwork3.pdf
diff --git a/usth/MATH2.4/labwork/3/subst.m b/usth/MATH2.4/labwork/3/subst.m
new file mode 100644
index 0000000..5c6e783
--- /dev/null
+++ b/usth/MATH2.4/labwork/3/subst.m
@@ -0,0 +1,17 @@
+function x = subst (aug)
+  [m n] = size (aug);
+  x = zeros (m, 1);
+  if (istril (aug))
+    x(m) = aug(1, n) / aug(1, m);
+    for k = 2 : m
+      x(k) = (aug(k, n) - aug(k, k+1:m)*x(1:k-1)) / aug(k, k);
+    end
+  elseif (istriu (aug))
+    x(m) = aug(m, n) / aug(m, m);
+    for k = m - 1 : -1 : 1
+      x(k) = (aug(k, n) - aug(k, k+1:m)*x(k+1:m)) / aug(k, k);
+    end
+  else
+    error ('aug must be a triangular matrix');
+  end
+end