https://mooseframework.inl.gov
BrentsMethod.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "BrentsMethod.h"
11 #include "MooseError.h"
12 #include "Conversion.h"
13 
14 namespace BrentsMethod
15 {
16 void
17 bracket(std::function<Real(Real)> const & f, Real & x1, Real & x2)
18 {
19  Real f1, f2;
20  // Factor to scale the interval by
21  Real factor = 1.6;
22  // Maximum number of times to attempt bracketing the interval
23  unsigned int n = 50;
24  // Small positive value to keep guess above
25  Real eps = 1.0e-10;
26 
27  // If the initial guesses are identical
28  if (x1 == x2)
29  throw MooseException("Bad initial range (0) used in BrentsMethod::bracket");
30 
31  f1 = f(x1);
32  f2 = f(x2);
33 
34  if (f1 * f2 > 0.0)
35  {
36  unsigned int iter = 0;
37  std::stringstream debug_ss;
38  while (f1 * f2 > 0.0)
39  {
40 #ifdef DEBUG
41  debug_ss << " iteration " << iter << ": (x1,x2) = (" << x1 << "," << x2 << "), (f1,f2) = ("
42  << f1 << "," << f2 << ")\n";
43 #endif
44  if (std::abs(f1) < std::abs(f2))
45  {
46  x1 += factor * (x1 - x2);
47  x1 = (x1 < eps ? eps : x1);
48  f1 = f(x1);
49  }
50  else
51  {
52  x2 += factor * (x2 - x1);
53  x2 = (x2 < eps ? eps : x2);
54  f2 = f(x2);
55  }
57  iter++;
58  if (iter >= n)
59  throw MooseException("No bracketing interval found by BrentsMethod::bracket after " +
60  Moose::stringify(n) + " iterations.\n" + debug_ss.str());
61  }
62  }
63 }
64 
65 Real
66 root(std::function<Real(Real)> const & f, Real x1, Real x2, Real tol)
67 {
68  Real a = x1, b = x2, c = x2, d = 0.0, e = 0.0, min1, min2;
69  Real fa = f(a);
70  Real fb = f(b);
71  Real fc, p, q, r, s, tol1, xm = 0;
72  unsigned int iter_max = 100;
73  Real eps = 1.0e-12;
74 
75  if (fa * fb > 0.0)
76  throw MooseException("Root must be bracketed in BrentsMethod::root");
77 
78  fc = fb;
79  std::stringstream debug_ss;
80  for (unsigned int i = 1; i <= iter_max; ++i)
81  {
82 #ifdef DEBUG
83  debug_ss << " iteration " << i << ": dx = " << xm << ", x = " << b << ", f(x) = " << fb
84  << "\n";
85 #endif
86  if (fb * fc > 0.0)
87  {
88  // Rename a,b and c and adjust bounding interval d
89  c = a;
90  fc = fa;
91  d = b - a;
92  e = d;
93  }
94  if (std::abs(fc) < std::abs(fb))
95  {
96  a = b;
97  b = c;
98  c = a;
99  fa = fb;
100  fb = fc;
101  fc = fa;
102  }
103  // Convergence check tolerance
104  tol1 = 2.0 * eps * std::abs(b) + 0.5 * tol;
105  xm = 0.5 * (c - b);
106 
107  if (std::abs(xm) <= tol1 || fb == 0.0)
108  return b;
109 
110  if (std::abs(e) >= tol1 && std::abs(fa) > std::abs(fb))
111  {
112  // Attempt inverse quadratic interpolation
113  s = fb / fa;
114  if (a == c)
115  {
116  p = 2.0 * xm * s;
117  q = 1.0 - s;
118  }
119  else
120  {
121  q = fa / fc;
122  r = fb / fc;
123  p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
124  q = (q - 1.0) * (r - 1.0) * (s - 1.0);
125  }
126  // Check whether in bounds
127  if (p > 0.0)
128  q = -q;
129  p = std::abs(p);
130  min1 = 3.0 * xm * q - std::abs(tol1 * q);
131  min2 = std::abs(e * q);
132 
133  if (2.0 * p < (min1 < min2 ? min1 : min2))
134  {
135  // Accept interpolation
136  e = d;
137  d = p / q;
138  }
139  else
140  {
141  // Interpolation failed, use bisection
142  d = xm;
143  e = d;
144  }
145  }
146  else
147  {
148  // Bounds decreasing too slowly, use bisection
149  d = xm;
150  e = d;
151  }
152  // Move last best guess to a
153  a = b;
154  fa = fb;
155  // Evaluate new trial root
156  if (std::abs(d) > tol1)
157  b += d;
158  else
159  {
160  Real sgn = (xm >= 0.0 ? std::abs(tol1) : -std::abs(tol1));
161  b += sgn;
162  }
163 
164  fb = f(b);
165  }
166 
167  throw MooseException("Maximum number of iterations exceeded in BrentsMethod::root.\n" +
168  debug_ss.str());
169  return 0.0; // Should never get here
170 }
171 } // namespace BrentsMethod
const Real eps
const double tol
int sgn(T val)
The sign function.
Definition: Numerics.h:41
Real f(Real x)
Test function for Brents method.
Real root(std::function< Real(Real)> const &f, Real x1, Real x2, Real tol=1.0e-12)
Finds the root of a function using Brent&#39;s method.
Definition: BrentsMethod.C:66
std::string stringify(const T &t)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Brent&#39;s method is used to find the root of a function f(x), i.e., find x such that f(x) = 0...
Definition: BrentsMethod.h:26
void bracket(std::function< Real(Real)> const &f, Real &x1, Real &x2)
Function to bracket a root of a given function.
Definition: BrentsMethod.C:17