BracketRoot.java

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.rootfinding;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.function.Function1D;

/**
 * Class that brackets single root of a function. For a 1-D function ({@link com.opengamma.analytics.math.function.Function1D}) $f(x)$,
 * initial values for the interval, $x_1$ and $x_2$, are supplied.
 * <p>
 * A root is assumed to be bracketed if $f(x_1)f(x_2) < 0$. If this condition is not satisfied, then either
 * $|f(x_1)| < |f(x_2)|$, in which case the lower value $x_1$ is shifted in the negative $x$ direction, or
 * the upper value $x_2$ is shifted in the positive $x$ direction. The amount by which to shift is the difference between
 * the two $x$ values multiplied by a constant ratio (1.6). If a root is not bracketed after 50 attempts, an exception is thrown.
 */
public class BracketRoot {
  private static final double RATIO = 1.6;
  private static final int MAX_STEPS = 50;

  /**
   * @param f The function, not null
   * @param xLower Initial value of lower bracket
   * @param xUpper Initial value of upper bracket
   * @return The bracketed points as an array, where the first element is the lower bracket and the second the upper bracket.
   * @throws MathException If a root is not bracketed in 50 attempts.
   */
  public double[] getBracketedPoints(final Function1D<Double, Double> f, final double xLower, final double xUpper) {
    Validate.notNull(f, "f");
    double x1 = xLower;
    double x2 = xUpper;
    double f1 = 0;
    double f2 = 0;
    f1 = f.evaluate(x1);
    f2 = f.evaluate(x2);
    if (Double.isNaN(f1)) {
      throw new MathException("Failed to bracket root: function invalid at x = " + x1 + " f(x) = " + f1);
    }
    if (Double.isNaN(f2)) {
      throw new MathException("Failed to bracket root: function invalid at x = " + x2 + " f(x) = " + f2);
    }

    for (int count = 0; count < MAX_STEPS; count++) {
      if (f1 * f2 < 0) {
        return new double[] {x1, x2 };
      }
      if (Math.abs(f1) < Math.abs(f2)) {
        x1 += RATIO * (x1 - x2);
        f1 = f.evaluate(x1);
        if (Double.isNaN(f1)) {
          throw new MathException("Failed to bracket root: function invalid at x = " + x1 + " f(x) = " + f1);
        }
      } else {
        x2 += RATIO * (x2 - x1);
        f2 = f.evaluate(x2);
        if (Double.isNaN(f2)) {
          throw new MathException("Failed to bracket root: function invalid at x = " + x2 + " f(x) = " + f2);
        }
      }
    }
    throw new MathException("Failed to bracket root");
  }

  public double[] getBracketedPoints(final Function1D<Double, Double> f, final double xLower, final double xUpper, final double minX, final double maxX) {
    Validate.notNull(f, "f");
    Validate.isTrue(xLower >= minX, "xLower < minX");
    Validate.isTrue(xUpper <= maxX, "xUpper < maxX");
    double x1 = xLower;
    double x2 = xUpper;
    double f1 = 0;
    double f2 = 0;
    boolean lowerLimitReached = false;
    boolean upperLimitReached = false;
    f1 = f.evaluate(x1);
    f2 = f.evaluate(x2);
    if (Double.isNaN(f1)) {
      throw new MathException("Failed to bracket root: function invalid at x = " + x1 + " f(x) = " + f1);
    }
    if (Double.isNaN(f2)) {
      throw new MathException("Failed to bracket root: function invalid at x = " + x2 + " f(x) = " + f2);
    }
    for (int count = 0; count < MAX_STEPS; count++) {
      if (f1 * f2 <= 0) {
        return new double[] {x1, x2 };
      }
      if (lowerLimitReached && upperLimitReached) {
        throw new MathException("Failed to bracket root: no root found between minX and maxX");
      }
      if (Math.abs(f1) < Math.abs(f2) && !lowerLimitReached) {
        x1 += RATIO * (x1 - x2);
        if (x1 < minX) {
          x1 = minX;
          lowerLimitReached = true;
        }
        f1 = f.evaluate(x1);
        if (Double.isNaN(f1)) {
          throw new MathException("Failed to bracket root: function invalid at x = " + x1 + " f(x) = " + f1);
        }
      } else {
        x2 += RATIO * (x2 - x1);
        if (x2 > maxX) {
          x2 = maxX;
          upperLimitReached = true;
        }
        f2 = f.evaluate(x2);
        if (Double.isNaN(f2)) {
          throw new MathException("Failed to bracket root: function invalid at x = " + x2 + " f(x) = " + f2);
        }
      }
    }
    throw new MathException("Failed to bracket root: max iterations");
  }

}