com.opengamma.analytics.math.rootfinding.BracketRoot.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.rootfinding.BracketRoot.java

Source

/**
 * 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");
    }

}