document + improvize root finding algorithms

computer_oriented_statistical_methods/false_position.cpp

@@ -1,19 +1,53 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [false position
+ * method](https://en.wikipedia.org/wiki/Regula_falsi), also known as the Secant
+ * method
+ *
+ * Given two points \f$a\f$ and \f$b\f$ such that \f$f(a)<0\f$ and
+ * \f$f(b)>0\f$, then the \f$(i+1)^\text{th}\f$ approximation is given by: \f[
+ * x_{i+1} = \frac{a_i\cdot f(b_i) - b_i\cdot f(a_i)}{f(b_i) - f(a_i)}
+ * \f]
+ * For the next iteration, the interval is selected
+ * as: \f$[a,x]\f$ if \f$x>0\f$ or \f$[x,b]\f$ if \f$x<0\f$. The Process is
+ * continued till a close enough approximation is achieved.
+ *
+ * \see newton_raphson_method.cpp, bisection_method.cpp
+ */
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
+#include <limits>
 
-static float eq(float i) {
-    return (pow(i, 3) - (4 * i) - 9);  // origial equation
+#define EPSILON \
+    1e-6  // std::numeric_limits<double>::epsilon()  ///< system accuracy limit
+#define MAX_ITERATIONS 50000  ///< Maximum number of iterations to check
+
+/** define \f$f(x)\f$ to find root for
+ */
+static double eq(double i) {
+    return (std::pow(i, 3) - (4 * i) - 9);  // origial equation
 }
 
+/** get the sign of any given number */
+template <typename T>
+int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
+/** main function */
 int main() {
-    float a, b, z, c, m, n;
-    system("clear");
-    for (int i = 0; i < 100; i++) {
-        z = eq(i);
-        if (z >= 0) {
-            b = i;
-            a = --i;
+    double a = -1, b = 1, x, z, m, n, c;
+    int i;
+
+    // loop to find initial intervals a, b
+    for (int i = 0; i < MAX_ITERATIONS; i++) {
+        z = eq(a);
+        x = eq(b);
+        if (sgn(z) == sgn(x)) {  // same signs, increase interval
+            b++;
+            a--;
+        } else {  // if opposite signs, we got our interval
             break;
         }
     }
@@ -21,18 +55,20 @@ int main() {
     std::cout << "\nFirst initial: " << a;
     std::cout << "\nSecond initial: " << b;
 
-    for (int i = 0; i < 100; i++) {
-        float h, d;
+    for (i = 0; i < MAX_ITERATIONS; i++) {
         m = eq(a);
         n = eq(b);
+
         c = ((a * n) - (b * m)) / (n - m);
+
         a = c;
         z = eq(c);
-        if (z > 0 && z < 0.09) {  // stoping criteria
+
+        if (std::abs(z) < EPSILON) {  // stoping criteria
             break;
         }
     }
 
-    std::cout << "\n\nRoot: " << c;
+    std::cout << "\n\nRoot: " << c << "\t\tSteps: " << i << std::endl;
     return 0;
 }