BrentsMethod.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7
8 #include "BrentsMethod.h"
9 #include "MooseError.h"
10
11 namespace BrentsMethod
12 {
13 void
14 bracket(std::function<Real(Real)> const & f, Real & x1, Real & x2)
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 }
56
57 Real
58 root(std::function<Real(Real)> const & f, Real x1, Real x2, Real tol)
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 }
157 } // namespace BrentsMethod
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:58
static const double tol
Definition: XFEMFuncs.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:14