regula_falsi.h

Go to the documentation of this file.

#include <functional>
#include <stdexcept>
#include <iostream>
#include <math.h>
#include <string>
 
namespace chemmisol {
 
    template<typename T>
        class RegulaFalsi {
            private:
                T a_0;
                T b_0;
 
                std::function<T(const T&)> f;
 
                T solve_iter(T a, T b, T f_a, T f_b, std::size_t n) const;
                T solve_eps(T a, T b, T f_a, T f_b, T eps) const;
 
            public:
                RegulaFalsi(T a, T b, std::function<T(const T&)> f) :
                    a_0(a), b_0(b), f(f) {
                    }
 
                T solve_iter(std::size_t n) const;
 
                T solve_eps(T eps) const;
        };
 
    template<typename T>
        T RegulaFalsi<T>::solve_iter(std::size_t n) const {
            T f_a = f(a_0);
            T f_b = f(b_0);
            if(f_a*f_b > 0)
                throw std::logic_error(
                        " f(a_0) = " + std::to_string(f_a) +
                        " and f(b_0) " + std::to_string(f_b) +
                        " should be of opposite signs. (a_0=" +
                        std::to_string(a_0) + ", b_0=" + std::to_string(b_0) + ")"
                        );
            return solve_iter(a_0, b_0, f_a, f_b, n);
        }
 
    template<typename T>
        T RegulaFalsi<T>::solve_iter(T a, T b, T f_a, T f_b, std::size_t n) const {
            if(f_a == f_b)
                return a;
            T c = (a * f_b - b * f_a) / (f_b - f_a);
            if(n==0)
                return c;
            T f_c = f(c);
            if(f_c == 0.0)
                return c;
            if(f_c * f_a < 0)
                return solve_iter(a, c, f_a, f_c, n-1);
            return solve_iter(c, b, f_c, f_b, n-1);
        }
 
    template<typename T>
        T RegulaFalsi<T>::solve_eps(T eps) const {
            T f_a = f(a_0);
            T f_b = f(b_0);
            if(f_a*f_b > 0)
                throw std::logic_error(
                        "f(a_0) = " + std::to_string(f_a) +
                        " and f(b_0) " + std::to_string(f_b) +
                        " should be of opposite signs. (a_0=" +
                        std::to_string(a_0) + ", b_0=" + std::to_string(b_0) + ")"
                        );
            return solve_eps(a_0, b_0, f_a, f_b, eps);
        }
 
    template<typename T>
        T RegulaFalsi<T>::solve_eps(T a, T b, T f_a, T f_b, T eps) const {
            if(f_a == f_b)
                return a;
            T c = (a * f_b - b * f_a) / (f_b - f_a);
            T f_c = f(c);
            if(std::abs(f_c) < eps)
                return c;
            if(f_c * f_a < 0)
                return solve_eps(a, c, f_a, f_c, eps);
            return solve_eps(c, b, f_c, f_b, eps);
        }
}