brent_method_extrema.cpp File Reference

Find real extrema of a univariate real function in a given interval using Brent's method. More...

#include <cassert>
#include <cmath>
#include <functional>
#include <iostream>
#include <limits>

Include dependency graph for brent_method_extrema.cpp:

Macros
#define	_USE_MATH_DEFINES
	required for MS Visual C++
#define	EPSILON
	system accuracy limit

Functions
double	get_minima (const std::function< double(double)> &f, double lim_a, double lim_b)
	Get the real root of a function in the given interval.
void	test1 ()
	Test function to find root for the function \(f(x)= (x-2)^2\) in the interval \([1,5]\) Expected result = 2.
void	test2 ()
	Test function to find root for the function \(f(x)= x^{\frac{1}{x}}\) in the interval \([-2,10]\) Expected result: \(e\approx 2.71828182845904509\).
void	test3 ()
	Test function to find maxima for the function \(f(x)= \cos x\) in the interval \([0,12]\) Expected result: \(\pi\approx 3.14159265358979312\).
int	main ()

Detailed Description

Find real extrema of a univariate real function in a given interval using Brent's method.

Refer the algorithm discoverer's publication online and also associated book:

‍R. P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey, 1973

See also

golden_search_extrema.cpp

Author

Krishna Vedala

Macro Definition Documentation

EPSILON

#define EPSILON

Value:

    std::sqrt(  \
        std::numeric_limits<double>::epsilon())

system accuracy limit

#define EPSILON \
    std::sqrt(  \
        std::numeric_limits<double>::epsilon())  ///< system accuracy limit

Function Documentation

get_minima()

double get_minima	(	const std::function< double(double)> &	f,
		double	lim_a,
		double	lim_b )

Get the real root of a function in the given interval.

Parameters

f	function to get root for
lim_a	lower limit of search window
lim_b	upper limit of search window

Returns

root found in the interval

                                {
    uint32_t iters = 0;
 
    if (lim_a > lim_b) {
        std::swap(lim_a, lim_b);
    } else if (std::abs(lim_a - lim_b) <= EPSILON) {
        std::cerr << "Search range must be greater than " << EPSILON << "\n";
        return lim_a;
    }
 
    // golden ratio value
    const double M_GOLDEN_RATIO = (3.f - std::sqrt(5.f)) / 2.f;
 
    double v = lim_a + M_GOLDEN_RATIO * (lim_b - lim_a);
    double u, w = v, x = v;
    double fu, fv = f(v);
    double fw = fv, fx = fv;
 
    double mid_point = (lim_a + lim_b) / 2.f;
    double p = 0, q = 0, r = 0;
 
    double d, e = 0;
    double tolerance, tolerance2;
 
    do {
        mid_point = (lim_a + lim_b) / 2.f;
        tolerance = EPSILON * std::abs(x);
        tolerance2 = 2 * tolerance;
 
        if (std::abs(e) > tolerance2) {
            // fit parabola
            r = (x - w) * (fx - fv);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w) * r;
            q = 2.f * (q - r);
            if (q > 0)
                p = -p;
            else
                q = -q;
            r = e;
            e = d;
        }
 
        if (std::abs(p) < std::abs(0.5 * q * r) && p < q * (lim_b - x)) {
            // parabolic interpolation step
            d = p / q;
            u = x + d;
            if (u - lim_a < tolerance2 || lim_b - u < tolerance2)
                d = x < mid_point ? tolerance : -tolerance;
        } else {
            // golden section interpolation step
            e = (x < mid_point ? lim_b : lim_a) - x;
            d = M_GOLDEN_RATIO * e;
        }
 
        // evaluate not too close to x
        if (std::abs(d) >= tolerance)
            u = d;
        else if (d > 0)
            u = tolerance;
        else
            u = -tolerance;
        u += x;
        fu = f(u);
 
        // update variables
        if (fu <= fx) {
            if (u < x)
                lim_b = x;
            else
                lim_a = x;
            v = w;
            fv = fw;
            w = x;
            fw = fx;
            x = u;
            fx = fu;
        } else {
            if (u < x)
                lim_a = u;
            else
                lim_b = u;
            if (fu <= fw || x == w) {
                v = w;
                fv = fw;
                w = u;
                fw = fu;
            } else if (fu <= fv || v == x || v == w) {
                v = u;
                fv = fu;
            }
        }
 
        iters++;
    } while (std::abs(x - mid_point) > (tolerance - (lim_b - lim_a) / 2.f));
 
    std::cout << " (iters: " << iters << ") ";
 
    return x;
}

Here is the call graph for this function:

main()

int main

(

void

)

Main function

           {
    std::cout.precision(18);
 
    std::cout << "Computations performed with machine epsilon: " << EPSILON
              << "\n";
 
    test1();
    test2();
    test3();
 
    return 0;
}

Here is the call graph for this function:

test1()

void test1

(

)

Test function to find root for the function \(f(x)= (x-2)^2\) in the interval \([1,5]\)
Expected result = 2.

             {
    // define the function to minimize as a lambda function
    std::function<double(double)> f1 = [](double x) {
        return (x - 2) * (x - 2);
    };
 
    std::cout << "Test 1.... ";
 
    double minima = get_minima(f1, -1, 5);
 
    std::cout << minima << "...";
 
    assert(std::abs(minima - 2) < EPSILON);
    std::cout << "passed\n";
}

Here is the call graph for this function:

test2()

void test2

(

)

Test function to find root for the function \(f(x)= x^{\frac{1}{x}}\) in the interval \([-2,10]\)
Expected result: \(e\approx 2.71828182845904509\).

             {
    // define the function to maximize as a lambda function
    // since we are maximixing, we negated the function return value
    std::function<double(double)> func = [](double x) {
        return -std::pow(x, 1.f / x);
    };
 
    std::cout << "Test 2.... ";
 
    double minima = get_minima(func, -2, 5);
 
    std::cout << minima << " (" << M_E << ")...";
 
    assert(std::abs(minima - M_E) < EPSILON);
    std::cout << "passed\n";
}

Here is the call graph for this function:

test3()

void test3

(

)

Test function to find maxima for the function \(f(x)= \cos x\) in the interval \([0,12]\)
Expected result: \(\pi\approx 3.14159265358979312\).

             {
    // define the function to maximize as a lambda function
    // since we are maximixing, we negated the function return value
    std::function<double(double)> func = [](double x) { return std::cos(x); };
 
    std::cout << "Test 3.... ";
 
    double minima = get_minima(func, -4, 12);
 
    std::cout << minima << " (" << M_PI << ")...";
 
    assert(std::abs(minima - M_PI) < EPSILON);
    std::cout << "passed\n";
}

Here is the call graph for this function: