Dekker.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <functional>
#include <utility>
 
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Roots/Iteration.hpp"
 
namespace nc::roots
{
    //================================================================================
    // Class Description:
    class Dekker : public Iteration
    {
    public:
        //============================================================================
        // Method Description:
        Dekker(const double epsilon, std::function<double(double)> f) noexcept :
            Iteration(epsilon),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        Dekker(const double epsilon, const uint32 maxNumIterations, std::function<double(double)> f) noexcept :
            Iteration(epsilon, maxNumIterations),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        ~Dekker() override = default;
 
        //============================================================================
        // Method Description:
        double solve(double a, double b)
        {
            resetNumberOfIterations();
 
            double fa = f_(a);
            double fb = f_(b);
 
            checkAndFixAlgorithmCriteria(a, b, fa, fb);
 
            double lastB  = a;
            double lastFb = fa;
 
            while (std::fabs(fb) > epsilon_ && std::fabs(b - a) > epsilon_)
            {
                const double s = calculateSecant(b, fb, lastB, lastFb);
                const double m = calculateBisection(a, b);
 
                lastB = b;
 
                b = useSecantMethod(b, s, m) ? s : m;
 
                lastFb = fb;
                fb     = f_(b);
 
                if (fa * fb > 0 && fb * lastFb < 0)
                {
                    a = lastB;
                }
 
                fa = f_(a);
                checkAndFixAlgorithmCriteria(a, b, fa, fb);
 
                incrementNumberOfIterations();
            }
 
            return b;
        }
 
    private:
        //============================================================================
        const std::function<double(double)> f_;
 
        //============================================================================
        // Method Description:
        static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) noexcept
        {
            // Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
            if (std::fabs(fa) < std::fabs(fb))
            {
                std::swap(a, b);
                std::swap(fa, fb);
            }
        }
 
        //============================================================================
        // Method Description:
        static double calculateSecant(double b, double fb, double lastB, double lastFb) noexcept
        {
            // No need to check division by 0, in this case the method returns NAN which is taken care by
            // useSecantMethod method
            return b - fb * (b - lastB) / (fb - lastFb);
        }
 
        //============================================================================
        // Method Description:
        static double calculateBisection(double a, double b) noexcept
        {
            return 0.5 * (a + b);
        }
 
        //============================================================================
        // Method Description:
        static bool useSecantMethod(double b, double s, double m) noexcept
        {
            // Value s calculated by secant method has to be between m and b
            return (b > m && s > m && s < b) || (b < m && s > b && s < m);
        }
    };
} // namespace nc::roots