root_finders.h

Go to the documentation of this file.

#ifndef WXM_ROOT_FINDERS_H
#define WXM_ROOT_FINDERS_H
 
#include "warpxm/warpxm_config.h"
#include <cassert>
 
namespace wxm
{
namespace root_finders
{
 
// Finding root of function using bisection.
// fun(a) and fun(b)  must be of opposite sign
// This is a modified version of:
// //
// https://stackoverflow.com/questions/28267910/function-as-parameter-of-a-function-using-bisection-method-c
// instead of function pointer we have template
// template <class T>
// real bisection(const T &fun, real a, real b);
 
// template <class T>
// real bisection(const T &fun, real a, real b) {
//     static const real EPS = 1e-15; // 1×10^(-15)
//     real fa = fun(a), fb = fun(b);
 
//  std::cout << "initially a = " << a << " b = " << b << " f(a) = " << fa << " f(b) = "
// << fb << std::endl;
 
//     // if either f(a) or f(b) are the root, return that
//     // nothing else to do
//     if (fa == 0) return a;
//     if (fb == 0) return b;
 
//     // this method only works if the signs of f(a) and f(b)
//     // are different. so just assert that
//     assert(fa * fb < 0); // 8.- macro assert from header cassert.
 
//     do {
//         // calculate fun at the midpoint of a,b
//         // if that's the root, we're done
 
//         // this line is awful, never write code like this...
//         //if ((f = fun((s = (a + b) / 2))) == 0) break;
 
//         // prefer:
//         real midpt = (a + b) / 2;
//         real fmid = fun(midpt);
//      std::cout << "(midpt,fmid) = (" << midpt << ", " << fmid << ")" << std::endl;
 
//         if (fmid == 0) return midpt;
 
//         // adjust our bounds to either [a,midpt] or [midpt,b]
//         // based on where fmid ends up being. I'm pretty
//         // sure the code in the question is wrong, so I fixed it
//         if (fa * fmid < 0) { // fmid, not f1
//             fb = fmid;
//             b = midpt;
//         }
//         else {
//             fa = fmid;
//             a = midpt;
//         }
//     } while (fabs(b-a) > EPS); // only loop while
//                          // a and b are sufficiently far
//                          // apart
 
//  std::cout << "a = " << a << " b = " << b << std::endl;
//     return (a + b) / 2;  // approximation
// }
 
// This is a modified version of:
// // https://www.geeksforgeeks.org/program-for-bisection-method/
// instead of function pointer we have template
// also there is a check that b >= a for the interval and if not they are swapped
 
template<class T> real bisection(const T& fun, real a, real b)
{
    // make sure a < b
 
    if (a > b)
    {
        int temp = b;
        b = a;
        a = temp;
    }
    // std::cout << "initially a = " << a << " b = " << b << " f(a) = " << fun(a) << "
    // f(b) = " << fun(b) << std::endl;
    static const real epsilon = 1e-15; // 1×10^(-15)
    if (fun(a) * fun(b) >= 0)
    {
        // throw WxExcept("wxm::root_finders::bisection You have not assumed correct
        // bounds a (=%f) and b\n",a);
        WxExcept wxe("wxm::root_finders::bisection You have not assumed correct bounds");
 
        throw wxe << " a = " << a << " b = " << b << std::endl;
        // return;
    }
 
    // std::cout << "initially (b-a) = " << (b-a) << std::endl;
    real c = a;
    // while (fabs(b-a) >= epsilon)
    while (b - a >= epsilon)
    {
        // std::cout << "(b-a) = " << (b-a) << std::endl;
        // Find middle point
        c = (a + b) / 2;
 
        // Check if middle point is root
        if (fun(c) == 0.0)
            return c;
 
        // Decide the side to repeat the steps
        else if (fun(c) * fun(a) < 0)
            b = c;
        else
            a = c;
    }
    // std::cout << "The value of root is : " << c << " and f(c) = " << fun(c) <<
    // std::endl;
 
    return c;
}
 
} // namespace root_finders
} // namespace wxm
 
#endif // WXM_ROOT_FINDERS_H