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 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)
 create mode 100644 usth/MATH2.2/labwork/1/bisect.m

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

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
-- 
cgit 1.4.1