From e461df7573c2b7b7e26c965d8cf2d8e175d67378 Mon Sep 17 00:00:00 2001
From: Nguyễn Gia Phong <vn.mcsinyx@gmail.com>
Date: Mon, 16 Dec 2019 21:31:18 +0700
Subject: [usth/MATH2.2] Numerical Methods

The future starts now.
---
 usth/MATH2.2/labwork/1/bisect.m     |  36 +++++++++++++++++++++
 usth/MATH2.2/labwork/1/intro        |  62 ++++++++++++++++++++++++++++++++++++
 usth/MATH2.2/labwork/1/labwork1.pdf | Bin 0 -> 165522 bytes
 usth/MATH2.2/labwork/1/ratpoison.m  |  15 +++++++++
 usth/MATH2.2/labwork/1/ratpoison.py |   8 +++++
 5 files changed, 121 insertions(+)
 create mode 100644 usth/MATH2.2/labwork/1/bisect.m
 create mode 100755 usth/MATH2.2/labwork/1/intro
 create mode 100644 usth/MATH2.2/labwork/1/labwork1.pdf
 create mode 100644 usth/MATH2.2/labwork/1/ratpoison.m
 create mode 100644 usth/MATH2.2/labwork/1/ratpoison.py

(limited to 'usth/MATH2.2/labwork/1')

diff --git a/usth/MATH2.2/labwork/1/bisect.m b/usth/MATH2.2/labwork/1/bisect.m
new file mode 100644
index 0000000..1d9e99f
--- /dev/null
+++ b/usth/MATH2.2/labwork/1/bisect.m
@@ -0,0 +1,36 @@
+function [x fx ea iter] = bisect (f, xl, xu, es = 0.0001, maxit = 50)
+  % uses bisection method to find the root of f,
+  % with xl and xu being lower and upper guesses,
+  % es being desired relative error
+  % and maxit being maximum allowable iterations
+  % to return the real root x, fx = f(x),
+  % approximate relative error ea (%)
+  % and number of iterations iter
+  nargin < 3 && error ('bisect requires at least 3 arguments');
+  [fl fu iter] = deal (f (xl), f (xu), 0);
+  if (fl == 0)
+    [x fx ea] = deal (xl, 0, 0);
+    return;
+  elseif (fu == 0)
+    [x fx ea] = deal (xu, 0, 0);
+    return;
+  end
+  fl * fu < 0 || error ('no sign change');
+
+  [x ea] = deal (xl, 100);
+  while (ea > es && iter++ < maxit) % yes, I use Octave only
+    [xold x] = deal (x, (xl + xu) / 2);
+    fx = f (x);
+    if (fx == 0)
+      ea = 0;
+      break;
+    elseif (x)
+      ea = abs ((x - xold) / x) * 100;
+    end
+    if (f (xl) * fx < 0)
+      xu = x;
+    else
+      xl = x;
+    end
+  end
+end
diff --git a/usth/MATH2.2/labwork/1/intro b/usth/MATH2.2/labwork/1/intro
new file mode 100755
index 0000000..555f874
--- /dev/null
+++ b/usth/MATH2.2/labwork/1/intro
@@ -0,0 +1,62 @@
+#!/usr/bin/env octave
+% Exercise I.1.1
+x = 1 : 5
+printf ("x (3) + x (2) = %f\n", x (3) + x (2));
+
+% Exercise I.1.2
+x = x'
+
+% Exercise I.2.1
+A = zeros (3, 3)
+A += eye (3)
+A += [0 0 0; 2 0 0; 3 4 0]
+B = A'
+A + B
+A - B
+A * B  % matmul
+A .* B % element-wise mul
+A(1:2, 1:2)
+A(1:3)
+A(1:3, 1) % col
+A(:, 1) % col
+
+% Exercise I.3
+function avg = mean (v)
+  avg = sum (v) / numel (v);
+end
+mean (x)
+
+% Exercise I.4
+function price = ebill (kwh)
+  l = [50 50 100 100 100 Inf];
+  p = [1.484 1.533 1.786 2.242 2.503 2.587];
+  for i = 1 : 6
+    if (kwh > l(i))
+      kwh -= l(i);
+    else
+      l(i) = kwh;
+      kwh = 0;
+    end
+  end
+  price = sum (l .* p);
+end
+
+% Exercise I.5.1
+a = (@(x) x*x') (1 : 1000)
+b = sum (arrayfun (@(x) (-1)^x / (1 + x*2), [0 : 501]))
+c = sum (arrayfun (@(x) 1 / (1 + x*2)^2 / (3 + x*2)^2, [0 : 420])) - (pi^2 - 8)/16
+
+% Exercise I.5.2
+function p = permutations (n, k)
+  p = factorial (n) / factorial (n - k);
+end
+function c = combinations (n, k)
+  c = permutations (n, k) / factorial (k);
+end
+
+% Exercise I.6
+log (2) / log (1.1)
+
+% Exercise I.7
+x = linspace(-5, 5);
+plot (x, x + 1, 'color', 'green');
diff --git a/usth/MATH2.2/labwork/1/labwork1.pdf b/usth/MATH2.2/labwork/1/labwork1.pdf
new file mode 100644
index 0000000..d8ce2d9
Binary files /dev/null and b/usth/MATH2.2/labwork/1/labwork1.pdf differ
diff --git a/usth/MATH2.2/labwork/1/ratpoison.m b/usth/MATH2.2/labwork/1/ratpoison.m
new file mode 100644
index 0000000..f2db823
--- /dev/null
+++ b/usth/MATH2.2/labwork/1/ratpoison.m
@@ -0,0 +1,15 @@
+function [x fx ea i] = ratpoison (f, df, x0, es = 0.00000001, imax = 20)
+  nargin < 2 && error ('ratpoison requires at least 2 ingredients');
+  [x fx dfx ea i] = deal (x0, f (x0), df (x0), 1, 0);
+  while (ea > es && i++ < imax)
+    [xold x] = deal (x, x - fx/dfx);
+    [fx dfx] = deal (f (x), df (x));
+    if (fx == 0)
+      ea = 0;
+      break;
+    elseif (x)
+      % just drop the percent BS
+      ea = abs ((x - xold) / x);
+    end
+  end
+end
diff --git a/usth/MATH2.2/labwork/1/ratpoison.py b/usth/MATH2.2/labwork/1/ratpoison.py
new file mode 100644
index 0000000..16aa65d
--- /dev/null
+++ b/usth/MATH2.2/labwork/1/ratpoison.py
@@ -0,0 +1,8 @@
+def ratpoison(f, df, x, es=10**-8, imax=20):
+    ea, i = 1, 0
+    while ea > es and i < imax:
+        i += 1
+        xold, x = x, x - f(x)/df(x)
+        if f(x) == 0: return x, 0, 0, i
+        if x: ea = abs((x - xold) / x)
+    return x, f(x), ea, i
-- 
cgit 1.4.1