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

◆ bracket()

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 17 of file BrentsMethod.C.

Referenced by IdealRealGasMixtureFluidProperties::p_T_from_v_e(), CaloricallyImperfectGas::rho_from_p_s(), HelmholtzFluidProperties::rho_from_p_T(), CO2FluidProperties::rho_from_p_T(), CaloricallyImperfectGas::setupLookupTables(), GeochemicalModelInterrogator::solveForT(), IdealRealGasMixtureFluidProperties::T_from_p_v(), TEST(), and IdealRealGasMixtureFluidProperties::v_from_p_T().

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  while (f1 * f2 > 0.0)
38  {
39  if (std::abs(f1) < std::abs(f2))
40  {
41  x1 += factor * (x1 - x2);
42  x1 = (x1 < eps ? eps : x1);
43  f1 = f(x1);
44  }
45  else
46  {
47  x2 += factor * (x2 - x1);
48  x2 = (x2 < eps ? eps : x2);
49  f2 = f(x2);
50  }
52  iter++;
53  if (iter >= n)
54  throw MooseException("No bracketing interval found by BrentsMethod::bracket after " +
55  Moose::stringify(n) + " iterations");
56  }
57  }
58 }
const Real eps
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Real f(Real x)
Test function for Brents method.
std::string stringify(const T &t)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ root()

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 61 of file BrentsMethod.C.

Referenced by ParameterStudyAction::inferMultiAppMode(), IdealRealGasMixtureFluidProperties::p_T_from_v_e(), PolycrystalHex::precomputeGrainStructure(), CaloricallyImperfectGas::rho_from_p_s(), HelmholtzFluidProperties::rho_from_p_T(), CO2FluidProperties::rho_from_p_T(), CaloricallyImperfectGas::setupLookupTables(), IdealRealGasMixtureFluidProperties::T_from_p_v(), TEST(), testExceptionMessage(), and IdealRealGasMixtureFluidProperties::v_from_p_T().

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