feat: add bisection method

numerical_methods/bisection_method.c

@@ -0,0 +1,106 @@
+/**
+ * @file
+ * @brief In mathematics, the [Bisection
+ * Method](https://en.wikipedia.org/wiki/Bisection_method) is a root-finding
+ * method that applies to any continuous function for which one knows two values
+ * with opposite signs.
+ * @details
+ * The method consists of repeatedly bisecting the interval
+ * defined by the two values and then selecting the subinterval in which the
+ * function changes sign, and therefore must contain a root. It is a very
+ * simple and robust method, but it is also relatively slow. Because of this,
+ * it is often used to obtain a rough approximation to a solution which is
+ * then used as a starting point for more rapidly converging methods.
+ * @author [Aybars Nazlica](https://github.com/aybarsnazlica)
+ */
+
+#include <assert.h>  /// for assert
+#include <math.h>    /// for fabs
+
+#define EPSILON 0.0001  // a small positive infinitesimal quantity
+#define NMAX 50         // maximum number of iterations
+
+/**
+ * @brief Function to check if two values have the same sign
+ * @param a First value
+ * @param b Second value
+ * @returns `sign`
+ */
+double sign(double a, double b)
+{
+    return (a > 0 && b > 0) + (a < 0 && b < 0) - (a > 0 && b < 0) -
+           (a < 0 && b > 0);
+}
+
+/**
+ * @brief Continuous function for root-finding
+ * @param x Real input variable
+ * @returns `function evaluation result`
+ */
+double func(double x)
+{
+    return x * x * x + 2.0 * x - 10.0;  // f(x) = x**3 + 2x - 10
+}
+
+/**
+ * @brief Root-finding method for a continuous function given two values with
+ * opposite signs
+ * @param x_left Lower endpoint value of the interval
+ * @param x_right Upper endpoint value of the interval
+ * @param tolerance Error threshold
+ * @returns `root of the function` if bisection method succeed within the
+ * maximum number of iterations
+ * @returns `-1` if bisection method fails
+ */
+double bisection(double x_left, double x_right, double tolerance)
+{
+    int n = 1;      // step counter
+    double middle;  // midpoint
+
+    while (n <= NMAX)
+    {
+        middle = (x_left + x_right) / 2;  // bisect the interval
+        double error = middle - x_left;
+
+        if (fabs(func(middle)) < EPSILON || error < tolerance)
+        {
+            return middle;
+        }
+
+        if (sign(func(middle), func(x_left)) > 0.0)
+        {
+            x_left = middle;  // new lower endpoint
+        }
+        else
+        {
+            x_right = middle;  // new upper endpoint
+        }
+
+        n++;  // increase step counter
+    }
+    return -1;  // method failed (maximum number of steps exceeded)
+}
+
+/**
+ * @brief Self-test implementations
+ * @returns void
+ */
+static void test()
+{
+    /* Compares root value that is found by the bisection method within a given
+     * floating point error*/
+    assert(fabs(bisection(1.0, 2.0, 0.0001) - 1.847473) <
+           EPSILON);  // the algorithm works as expected
+    assert(fabs(bisection(100.0, 250.0, 0.0001) - 249.999928) <
+           EPSILON);  // the algorithm works as expected
+}
+
+/**
+ * @brief Main function
+ * @returns 0 on exit
+ */
+int main()
+{
+    test();  // run self-test implementations
+    return 0;
+}