Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
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 More...
 

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. More...
 
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:

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
ffunction to get root for
lim_alower limit of search window
lim_bupper limit of search window
Returns
root found in the interval
36  {
37  uint32_t iters = 0;
38 
39  if (lim_a > lim_b) {
40  std::swap(lim_a, lim_b);
41  } else if (std::abs(lim_a - lim_b) <= EPSILON) {
42  std::cerr << "Search range must be greater than " << EPSILON << "\n";
43  return lim_a;
44  }
45 
46  // golden ratio value
47  const double M_GOLDEN_RATIO = (3.f - std::sqrt(5.f)) / 2.f;
48 
49  double v = lim_a + M_GOLDEN_RATIO * (lim_b - lim_a);
50  double u, w = v, x = v;
51  double fu, fv = f(v);
52  double fw = fv, fx = fv;
53 
54  double mid_point = (lim_a + lim_b) / 2.f;
55  double p = 0, q = 0, r = 0;
56 
57  double d, e = 0;
58  double tolerance, tolerance2;
59 
60  do {
61  mid_point = (lim_a + lim_b) / 2.f;
62  tolerance = EPSILON * std::abs(x);
63  tolerance2 = 2 * tolerance;
64 
65  if (std::abs(e) > tolerance2) {
66  // fit parabola
67  r = (x - w) * (fx - fv);
68  q = (x - v) * (fx - fw);
69  p = (x - v) * q - (x - w) * r;
70  q = 2.f * (q - r);
71  if (q > 0)
72  p = -p;
73  else
74  q = -q;
75  r = e;
76  e = d;
77  }
78 
79  if (std::abs(p) < std::abs(0.5 * q * r) && p < q * (lim_b - x)) {
80  // parabolic interpolation step
81  d = p / q;
82  u = x + d;
83  if (u - lim_a < tolerance2 || lim_b - u < tolerance2)
84  d = x < mid_point ? tolerance : -tolerance;
85  } else {
86  // golden section interpolation step
87  e = (x < mid_point ? lim_b : lim_a) - x;
88  d = M_GOLDEN_RATIO * e;
89  }
90 
91  // evaluate not too close to x
92  if (std::abs(d) >= tolerance)
93  u = d;
94  else if (d > 0)
95  u = tolerance;
96  else
97  u = -tolerance;
98  u += x;
99  fu = f(u);
100 
101  // update variables
102  if (fu <= fx) {
103  if (u < x)
104  lim_b = x;
105  else
106  lim_a = x;
107  v = w;
108  fv = fw;
109  w = x;
110  fw = fx;
111  x = u;
112  fx = fu;
113  } else {
114  if (u < x)
115  lim_a = u;
116  else
117  lim_b = u;
118  if (fu <= fw || x == w) {
119  v = w;
120  fv = fw;
121  w = u;
122  fw = fu;
123  } else if (fu <= fv || v == x || v == w) {
124  v = u;
125  fv = fu;
126  }
127  }
128 
129  iters++;
130  } while (std::abs(x - mid_point) > (tolerance - (lim_b - lim_a) / 2.f));
131 
132  std::cout << " (iters: " << iters << ") ";
133 
134  return x;
135 }
Here is the call graph for this function:

◆ main()

int main ( void  )

Main function

204  {
205  std::cout.precision(18);
206 
207  std::cout << "Computations performed with machine epsilon: " << EPSILON
208  << "\n";
209 
210  test1();
211  test2();
212  test3();
213 
214  return 0;
215 }
Here is the call graph for this function:
test3
void test3()
Test function to find maxima for the function in the interval Expected result: .
Definition: brent_method_extrema.cpp:188
EPSILON
#define EPSILON
system accuracy limit
Definition: brent_method_extrema.cpp:23
std::sqrt
T sqrt(T... args)
std::cerr
std::swap
T swap(T... args)
test2
void test2()
Test function to find root for the function in the interval Expected result: .
Definition: brent_method_extrema.cpp:165
test1
void test1()
Test function to find root for the function in the interval Expected result = 2.
Definition: brent_method_extrema.cpp:143
std::numeric_limits