Secant.hpp

Go to the documentation of this file.

#pragma once
 
#include <cmath>
#include <functional>
#include <utility>
 
#include "NumCpp/Core/DtypeInfo.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Roots/Iteration.hpp"
 
namespace nc::roots
{
    //================================================================================
    // Class Description:
    class Secant : public Iteration
    {
    public:
        //============================================================================
        // Method Description:
        Secant(const double epsilon, std::function<double(double)> f) noexcept :
            Iteration(epsilon),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        Secant(const double epsilon, const uint32 maxNumIterations, std::function<double(double)> f) noexcept :
            Iteration(epsilon, maxNumIterations),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        ~Secant() override = default;
 
        //============================================================================
        // Method Description:
        double solve(double a, double b)
        {
            resetNumberOfIterations();
 
            if (f_(a) > f_(b))
            {
                std::swap(a, b);
            }
 
            double x      = b;
            double lastX  = a;
            double fx     = f_(b);
            double lastFx = f_(a);
 
            while (std::fabs(fx) >= epsilon_)
            {
                const double x_tmp = calculateX(x, lastX, fx, lastFx);
 
                lastFx = fx;
                lastX  = x;
                x      = x_tmp;
 
                fx = f_(x);
 
                incrementNumberOfIterations();
            }
 
            return x;
        }
 
    private:
        //============================================================================
        const std::function<double(double)> f_;
 
        //============================================================================
        // Method Description:
        static double calculateX(double x, double lastX, double fx, double lastFx) noexcept
        {
            const double functionDifference = fx - lastFx;
            return x - fx * (x - lastX) / functionDifference;
        }
    };
} // namespace nc::roots