BrentsMethod Namespace Reference

Brent's method is used to find the root of a function f(x), i.e., find x such that f(x) = 0. More...

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

Detailed Description

Brent's method is used to find the root of a function f(x), i.e., find x such that f(x) = 0.

First, brackets x1 and x2 are found such that f(x) changes sign between x1 and x2, implying that there is a root between the points.

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

	f	reference to function to find bracketing interval
[out]	x1	reference one bound
[out]	x2	reference to other bound

Increment counter

Increment counter

Definition at line 17 of file BrentsMethod.C.

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

{
  Real f1, f2;
  // Factor to scale the interval by
  Real factor = 1.6;
  // Maximum number of times to attempt bracketing the interval
  unsigned int n = 50;
  // Small positive value to keep guess above
  Real eps = 1.0e-10;

  // If the initial guesses are identical
  if (x1 == x2)
    throw MooseException("Bad initial range (0) used in BrentsMethod::bracket");

  f1 = f(x1);
  f2 = f(x2);

  if (f1 * f2 > 0.0)
  {
    unsigned int iter = 0;
    std::stringstream debug_ss;
    while (f1 * f2 > 0.0)
    {
#ifdef DEBUG
      debug_ss << "  iteration " << iter << ": (x1,x2) = (" << x1 << "," << x2 << "), (f1,f2) = ("
               << f1 << "," << f2 << ")\n";
#endif
      if (std::abs(f1) < std::abs(f2))
      {
        x1 += factor * (x1 - x2);
        x1 = (x1 < eps ? eps : x1);
        f1 = f(x1);
      }
      else
      {
        x2 += factor * (x2 - x1);
        x2 = (x2 < eps ? eps : x2);
        f2 = f(x2);
      }
      iter++;
      if (iter >= n)
        throw MooseException("No bracketing interval found by BrentsMethod::bracket after " +
                             Moose::stringify(n) + " iterations.\n" + debug_ss.str());
    }
  }
}

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

f	reference to function to find root of
x1	one end of bracketing interval
x2	other end of bracketing interval
tolerance	root finding tolerance (default is 1e-12)

Definition at line 66 of file BrentsMethod.C.

Referenced by InterWrapper1PhaseProblem::computeWijFromSolve(), SubChannel1PhaseProblem::computeWijFromSolve(), ParameterStudyAction::inferMultiAppMode(), InterWrapper1PhaseProblem::petscSnesSolver(), SubChannel1PhaseProblem::petscSnesSolver(), PolycrystalHex::precomputeGrainStructure(), CaloricallyImperfectGas::rho_from_p_s(), HelmholtzFluidProperties::rho_from_p_T(), CO2FluidProperties::rho_from_p_T(), CaloricallyImperfectGas::setupLookupTables(), TEST(), and testExceptionMessage().

{
  Real a = x1, b = x2, c = x2, d = 0.0, e = 0.0, min1, min2;
  Real fa = f(a);
  Real fb = f(b);
  Real fc, p, q, r, s, tol1, xm = 0;
  unsigned int iter_max = 100;
  Real eps = 1.0e-12;

  if (fa * fb > 0.0)
    throw MooseException("Root must be bracketed in BrentsMethod::root");

  fc = fb;
  std::stringstream debug_ss;
  for (unsigned int i = 1; i <= iter_max; ++i)
  {
#ifdef DEBUG
    debug_ss << "  iteration " << i << ": dx = " << xm << ", x = " << b << ", f(x) = " << fb
             << "\n";
#endif
    if (fb * fc > 0.0)
    {
      // Rename a,b and c and adjust bounding interval d
      c = a;
      fc = fa;
      d = b - a;
      e = d;
    }
    if (std::abs(fc) < std::abs(fb))
    {
      a = b;
      b = c;
      c = a;
      fa = fb;
      fb = fc;
      fc = fa;
    }
    // Convergence check tolerance
    tol1 = 2.0 * eps * std::abs(b) + 0.5 * tol;
    xm = 0.5 * (c - b);

    if (std::abs(xm) <= tol1 || fb == 0.0)
      return b;

    if (std::abs(e) >= tol1 && std::abs(fa) > std::abs(fb))
    {
      // Attempt inverse quadratic interpolation
      s = fb / fa;
      if (a == c)
      {
        p = 2.0 * xm * s;
        q = 1.0 - s;
      }
      else
      {
        q = fa / fc;
        r = fb / fc;
        p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
        q = (q - 1.0) * (r - 1.0) * (s - 1.0);
      }
      // Check whether in bounds
      if (p > 0.0)
        q = -q;
      p = std::abs(p);
      min1 = 3.0 * xm * q - std::abs(tol1 * q);
      min2 = std::abs(e * q);

      if (2.0 * p < (min1 < min2 ? min1 : min2))
      {
        // Accept interpolation
        e = d;
        d = p / q;
      }
      else
      {
        // Interpolation failed, use bisection
        d = xm;
        e = d;
      }
    }
    else
    {
      // Bounds decreasing too slowly, use bisection
      d = xm;
      e = d;
    }
    // Move last best guess to a
    a = b;
    fa = fb;
    // Evaluate new trial root
    if (std::abs(d) > tol1)
      b += d;
    else
    {
      Real sgn = (xm >= 0.0 ? std::abs(tol1) : -std::abs(tol1));
      b += sgn;
    }

    fb = f(b);
  }

  throw MooseException("Maximum number of iterations exceeded in BrentsMethod::root.\n" +
                       debug_ss.str());
  return 0.0; // Should never get here
}