www.mooseframework.org
Functions
BrentsMethod Namespace Reference

Functions

void bracket (std::function< Real(Real)> const &f, Real &x1, Real &x2)
 Function to bracket a root of a given function. More...
 
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's method. More...
 

Function Documentation

void BrentsMethod::bracket ( std::function< Real(Real)> const &  f,
Real &  x1,
Real &  x2 
)

Function to bracket a root of a given function.

Adapted from Numerical Recipes in C

Parameters
freference to function to find bracketing interval
[out]x1reference one bound
[out]x2reference to other bound

Increment counter

Increment counter

Definition at line 14 of file BrentsMethod.C.

Referenced by CO2FluidProperties::rho().

15 {
16  Real f1, f2;
17  // Factor to scale the interval by
18  Real factor = 1.6;
19  // Maximum number of times to attempt bracketing the interval
20  unsigned int n = 50;
21  // Small positive value to keep guess above
22  Real eps = 1.0e-10;
23 
24  // If the initial guesses are identical
25  if (x1 == x2)
26  mooseError("Bad initial range (0) used in BrentsMethod::bracket");
27 
28  f1 = f(x1);
29  f2 = f(x2);
30 
31  if (f1 * f2 > 0.0)
32  {
33  unsigned int iter = 0;
34  while (f1 * f2 > 0.0)
35  {
36  if (std::abs(f1) < std::abs(f2))
37  {
38  x1 += factor * (x1 - x2);
39  x1 = (x1 < eps ? eps : x1);
40  f1 = f(x1);
41  }
42  else
43  {
44  x2 += factor * (x2 - x1);
45  x2 = (x2 < eps ? eps : x2);
46  f2 = f(x2);
47  }
49  iter++;
50  if (iter >= n)
51  mooseError(
52  "No bracketing interval found by BrentsMethod::bracket after ", n, " iterations");
53  }
54  }
55 }
Real BrentsMethod::root ( std::function< Real(Real)> const &  f,
Real  x1,
Real  x2,
Real  tol = 1.0e-12 
)

Finds the root of a function using Brent's method.

Adapted from Numerical Recipes in C

Parameters
freference to function to find root of
x1one end of bracketing interval
x2other end of bracketing interval
toleranceroot finding tolerance (default is 1e-12)

Definition at line 58 of file BrentsMethod.C.

Referenced by PolycrystalHex::precomputeGrainStructure(), and CO2FluidProperties::rho().

59 {
60  Real a = x1, b = x2, c = x2, d = 0.0, e = 0.0, min1, min2;
61  Real fa = f(a);
62  Real fb = f(b);
63  Real fc, p, q, r, s, tol1, xm;
64  unsigned int iter_max = 100;
65  Real eps = 1.0e-12;
66 
67  if (fa * fb > 0.0)
68  mooseError("Root must be bracketed in BrentsMethod::root");
69 
70  fc = fb;
71  for (unsigned int i = 1; i <= iter_max; ++i)
72  {
73  if (fb * fc > 0.0)
74  {
75  // Rename a,b and c and adjust bounding interval d
76  c = a;
77  fc = fa;
78  d = b - a;
79  e = d;
80  }
81  if (std::abs(fc) < std::abs(fb))
82  {
83  a = b;
84  b = c;
85  c = a;
86  fa = fb;
87  fb = fc;
88  fc = fa;
89  }
90  // Convergence check tolerance
91  tol1 = 2.0 * eps * std::abs(b) + 0.5 * tol;
92  xm = 0.5 * (c - b);
93 
94  if (std::abs(xm) <= tol1 || fb == 0.0)
95  return b;
96 
97  if (std::abs(e) >= tol1 && std::abs(fa) > std::abs(fb))
98  {
99  // Attempt inverse quadratic interpolation
100  s = fb / fa;
101  if (a == c)
102  {
103  p = 2.0 * xm * s;
104  q = 1.0 - s;
105  }
106  else
107  {
108  q = fa / fc;
109  r = fb / fc;
110  p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
111  q = (q - 1.0) * (r - 1.0) * (s - 1.0);
112  }
113  // Check whether in bounds
114  if (p > 0.0)
115  q = -q;
116  p = std::abs(p);
117  min1 = 3.0 * xm * q - std::abs(tol1 * q);
118  min2 = std::abs(e * q);
119 
120  if (2.0 * p < (min1 < min2 ? min1 : min2))
121  {
122  // Accept interpolation
123  e = d;
124  d = p / q;
125  }
126  else
127  {
128  // Interpolation failed, use bisection
129  d = xm;
130  e = d;
131  }
132  }
133  else
134  {
135  // Bounds decreasing too slowly, use bisection
136  d = xm;
137  e = d;
138  }
139  // Move last best guess to a
140  a = b;
141  fa = fb;
142  // Evaluate new trial root
143  if (std::abs(d) > tol1)
144  b += d;
145  else
146  {
147  Real sgn = (xm >= 0.0 ? std::abs(tol1) : -std::abs(tol1));
148  b += sgn;
149  }
150 
151  fb = f(b);
152  }
153 
154  mooseError("Maximum number of iterations exceeded in BrentsMethod::root");
155  return 0.0; // Should never get here
156 }
static const double tol
Definition: XFEMFuncs.h:26