/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math3.analysis.solvers;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;

Utility routines for UnivariateSolver objects.
/** * Utility routines for {@link UnivariateSolver} objects. * */
public class UnivariateSolverUtils {
Class contains only static methods.
/** * Class contains only static methods. */
private UnivariateSolverUtils() {}
Convenience method to find a zero of a univariate real function. A default solver is used.
Params:
  • function – Function.
  • x0 – Lower bound for the interval.
  • x1 – Upper bound for the interval.
Throws:
Returns:a value where the function is zero.
/** * Convenience method to find a zero of a univariate real function. A default * solver is used. * * @param function Function. * @param x0 Lower bound for the interval. * @param x1 Upper bound for the interval. * @return a value where the function is zero. * @throws NoBracketingException if the function has the same sign at the * endpoints. * @throws NullArgumentException if {@code function} is {@code null}. */
public static double solve(UnivariateFunction function, double x0, double x1) throws NullArgumentException, NoBracketingException { if (function == null) { throw new NullArgumentException(LocalizedFormats.FUNCTION); } final UnivariateSolver solver = new BrentSolver(); return solver.solve(Integer.MAX_VALUE, function, x0, x1); }
Convenience method to find a zero of a univariate real function. A default solver is used.
Params:
  • function – Function.
  • x0 – Lower bound for the interval.
  • x1 – Upper bound for the interval.
  • absoluteAccuracy – Accuracy to be used by the solver.
Throws:
Returns:a value where the function is zero.
/** * Convenience method to find a zero of a univariate real function. A default * solver is used. * * @param function Function. * @param x0 Lower bound for the interval. * @param x1 Upper bound for the interval. * @param absoluteAccuracy Accuracy to be used by the solver. * @return a value where the function is zero. * @throws NoBracketingException if the function has the same sign at the * endpoints. * @throws NullArgumentException if {@code function} is {@code null}. */
public static double solve(UnivariateFunction function, double x0, double x1, double absoluteAccuracy) throws NullArgumentException, NoBracketingException { if (function == null) { throw new NullArgumentException(LocalizedFormats.FUNCTION); } final UnivariateSolver solver = new BrentSolver(absoluteAccuracy); return solver.solve(Integer.MAX_VALUE, function, x0, x1); }
Force a root found by a non-bracketing solver to lie on a specified side, as if the solver were a bracketing one.
Params:
  • maxEval – maximal number of new evaluations of the function (evaluations already done for finding the root should have already been subtracted from this number)
  • f – function to solve
  • bracketing – bracketing solver to use for shifting the root
  • baseRoot – original root found by a previous non-bracketing solver
  • min – minimal bound of the search interval
  • max – maximal bound of the search interval
  • allowedSolution – the kind of solutions that the root-finding algorithm may accept as solutions.
Throws:
Returns:a root approximation, on the specified side of the exact root
/** * Force a root found by a non-bracketing solver to lie on a specified side, * as if the solver were a bracketing one. * * @param maxEval maximal number of new evaluations of the function * (evaluations already done for finding the root should have already been subtracted * from this number) * @param f function to solve * @param bracketing bracketing solver to use for shifting the root * @param baseRoot original root found by a previous non-bracketing solver * @param min minimal bound of the search interval * @param max maximal bound of the search interval * @param allowedSolution the kind of solutions that the root-finding algorithm may * accept as solutions. * @return a root approximation, on the specified side of the exact root * @throws NoBracketingException if the function has the same sign at the * endpoints. */
public static double forceSide(final int maxEval, final UnivariateFunction f, final BracketedUnivariateSolver<UnivariateFunction> bracketing, final double baseRoot, final double min, final double max, final AllowedSolution allowedSolution) throws NoBracketingException { if (allowedSolution == AllowedSolution.ANY_SIDE) { // no further bracketing required return baseRoot; } // find a very small interval bracketing the root final double step = FastMath.max(bracketing.getAbsoluteAccuracy(), FastMath.abs(baseRoot * bracketing.getRelativeAccuracy())); double xLo = FastMath.max(min, baseRoot - step); double fLo = f.value(xLo); double xHi = FastMath.min(max, baseRoot + step); double fHi = f.value(xHi); int remainingEval = maxEval - 2; while (remainingEval > 0) { if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) { // compute the root on the selected side return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution); } // try increasing the interval boolean changeLo = false; boolean changeHi = false; if (fLo < fHi) { // increasing function if (fLo >= 0) { changeLo = true; } else { changeHi = true; } } else if (fLo > fHi) { // decreasing function if (fLo <= 0) { changeLo = true; } else { changeHi = true; } } else { // unknown variation changeLo = true; changeHi = true; } // update the lower bound if (changeLo) { xLo = FastMath.max(min, xLo - step); fLo = f.value(xLo); remainingEval--; } // update the higher bound if (changeHi) { xHi = FastMath.min(max, xHi + step); fHi = f.value(xHi); remainingEval--; } } throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING, xLo, xHi, fLo, fHi, maxEval - remainingEval, maxEval, baseRoot, min, max); }
This method simply calls bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations) with q and r set to 1.0 and maximumIterations set to Integer.MAX_VALUE.

Note: this method can take Integer.MAX_VALUE iterations to throw a ConvergenceException. Unless you are confident that there is a root between lowerBound and upperBound near initial, it is better to use bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations), explicitly specifying the maximum number of iterations.

Params:
  • function – Function.
  • initial – Initial midpoint of interval being expanded to bracket a root.
  • lowerBound – Lower bound (a is never lower than this value)
  • upperBound – Upper bound (b never is greater than this value).
Throws:
Returns:a two-element array holding a and b.
/** * This method simply calls {@link #bracket(UnivariateFunction, double, double, double, * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)} * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}. * <p> * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE} * iterations to throw a {@code ConvergenceException.} Unless you are * confident that there is a root between {@code lowerBound} and * {@code upperBound} near {@code initial}, it is better to use * {@link #bracket(UnivariateFunction, double, double, double, double,double, int) * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}, * explicitly specifying the maximum number of iterations.</p> * * @param function Function. * @param initial Initial midpoint of interval being expanded to * bracket a root. * @param lowerBound Lower bound (a is never lower than this value) * @param upperBound Upper bound (b never is greater than this * value). * @return a two-element array holding a and b. * @throws NoBracketingException if a root cannot be bracketted. * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}. * @throws NullArgumentException if {@code function} is {@code null}. */
public static double[] bracket(UnivariateFunction function, double initial, double lowerBound, double upperBound) throws NullArgumentException, NotStrictlyPositiveException, NoBracketingException { return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, Integer.MAX_VALUE); }
This method simply calls bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations) with q and r set to 1.0.
Params:
  • function – Function.
  • initial – Initial midpoint of interval being expanded to bracket a root.
  • lowerBound – Lower bound (a is never lower than this value).
  • upperBound – Upper bound (b never is greater than this value).
  • maximumIterations – Maximum number of iterations to perform
Throws:
Returns:a two element array holding a and b.
/** * This method simply calls {@link #bracket(UnivariateFunction, double, double, double, * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)} * with {@code q} and {@code r} set to 1.0. * @param function Function. * @param initial Initial midpoint of interval being expanded to * bracket a root. * @param lowerBound Lower bound (a is never lower than this value). * @param upperBound Upper bound (b never is greater than this * value). * @param maximumIterations Maximum number of iterations to perform * @return a two element array holding a and b. * @throws NoBracketingException if the algorithm fails to find a and b * satisfying the desired conditions. * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}. * @throws NullArgumentException if {@code function} is {@code null}. */
public static double[] bracket(UnivariateFunction function, double initial, double lowerBound, double upperBound, int maximumIterations) throws NullArgumentException, NotStrictlyPositiveException, NoBracketingException { return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations); }
This method attempts to find two values a and b satisfying
  • lowerBound <= a < initial < b <= upperBound
  • f(a) * f(b) <= 0
If f is continuous on [a,b], this means that a and b bracket a root of f.

The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing values of k, where \( l_k = max(lower, initial - \delta_k) \), \( u_k = min(upper, initial + \delta_k) \), using recurrence \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \). The algorithm stops when one of the following happens:

  • at least one positive and one negative value have been found -- success!
  • both endpoints have reached their respective limits -- NoBracketingException
  • maximumIterations iterations elapse -- NoBracketingException

If different signs are found at first iteration (k=1), then the returned interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later iteration k>1, then the returned interval will be either \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called with these parameters will therefore start with the smallest bracketing interval known at this step.

Interval expansion rate is tuned by changing the recurrence parameters r and q. When the multiplicative factor r is set to 1, the sequence is a simple arithmetic sequence with linear increase. When the multiplicative factor r is larger than 1, the sequence has an asymptotically exponential rate. Note than the additive parameter q should never be set to zero, otherwise the interval would degenerate to the single initial point for all values of k.

As a rule of thumb, when the location of the root is expected to be approximately known within some error margin, r should be set to 1 and q should be set to the order of magnitude of the error margin. When the location of the root is really a wild guess, then r should be set to a value larger than 1 (typically 2 to double the interval length at each iteration) and q should be set according to half the initial search interval length.

As an example, if we consider the trivial function f(x) = 1 - x and use initial = 4, r = 1, q = 2, the algorithm will compute f(4-2) = f(2) = -1 and f(4+2) = f(6) = -5 for k = 1, then f(4-4) = f(0) = +1 and f(4+4) = f(8) = -7 for k = 2. Then it will return the interval [0, 2] as the smallest one known to be bracketing the root. As shown by this example, the initial value (here 4) may lie outside of the returned bracketing interval.

Params:
  • function – function to check
  • initial – Initial midpoint of interval being expanded to bracket a root.
  • lowerBound – Lower bound (a is never lower than this value).
  • upperBound – Upper bound (b never is greater than this value).
  • q – additive offset used to compute bounds sequence (must be strictly positive)
  • r – multiplicative factor used to compute bounds sequence
  • maximumIterations – Maximum number of iterations to perform
Throws:
Returns:a two element array holding the bracketing values.
/** * This method attempts to find two values a and b satisfying <ul> * <li> {@code lowerBound <= a < initial < b <= upperBound} </li> * <li> {@code f(a) * f(b) <= 0} </li> * </ul> * If {@code f} is continuous on {@code [a,b]}, this means that {@code a} * and {@code b} bracket a root of {@code f}. * <p> * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing * values of k, where \( l_k = max(lower, initial - \delta_k) \), * \( u_k = min(upper, initial + \delta_k) \), using recurrence * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \). * The algorithm stops when one of the following happens: <ul> * <li> at least one positive and one negative value have been found -- success!</li> * <li> both endpoints have reached their respective limits -- NoBracketingException </li> * <li> {@code maximumIterations} iterations elapse -- NoBracketingException </li></ul> * <p> * If different signs are found at first iteration ({@code k=1}), then the returned * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later * iteration {@code k>1}, then the returned interval will be either * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called * with these parameters will therefore start with the smallest bracketing interval known * at this step. * </p> * <p> * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r} * is larger than 1, the sequence has an asymptotically exponential rate. Note than the * additive parameter {@code q} should never be set to zero, otherwise the interval would * degenerate to the single initial point for all values of {@code k}. * </p> * <p> * As a rule of thumb, when the location of the root is expected to be approximately known * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the * order of magnitude of the error margin. When the location of the root is really a wild guess, * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval * length at each iteration) and {@code q} should be set according to half the initial * search interval length. * </p> * <p> * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root. * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned * bracketing interval. * </p> * @param function function to check * @param initial Initial midpoint of interval being expanded to * bracket a root. * @param lowerBound Lower bound (a is never lower than this value). * @param upperBound Upper bound (b never is greater than this * value). * @param q additive offset used to compute bounds sequence (must be strictly positive) * @param r multiplicative factor used to compute bounds sequence * @param maximumIterations Maximum number of iterations to perform * @return a two element array holding the bracketing values. * @exception NoBracketingException if function cannot be bracketed in the search interval */
public static double[] bracket(final UnivariateFunction function, final double initial, final double lowerBound, final double upperBound, final double q, final double r, final int maximumIterations) throws NoBracketingException { if (function == null) { throw new NullArgumentException(LocalizedFormats.FUNCTION); } if (q <= 0) { throw new NotStrictlyPositiveException(q); } if (maximumIterations <= 0) { throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations); } verifySequence(lowerBound, initial, upperBound); // initialize the recurrence double a = initial; double b = initial; double fa = Double.NaN; double fb = Double.NaN; double delta = 0; for (int numIterations = 0; (numIterations < maximumIterations) && (a > lowerBound || b < upperBound); ++numIterations) { final double previousA = a; final double previousFa = fa; final double previousB = b; final double previousFb = fb; delta = r * delta + q; a = FastMath.max(initial - delta, lowerBound); b = FastMath.min(initial + delta, upperBound); fa = function.value(a); fb = function.value(b); if (numIterations == 0) { // at first iteration, we don't have a previous interval // we simply compare both sides of the initial interval if (fa * fb <= 0) { // the first interval already brackets a root return new double[] { a, b }; } } else { // we have a previous interval with constant sign and expand it, // we expect sign changes to occur at boundaries if (fa * previousFa <= 0) { // sign change detected at near lower bound return new double[] { a, previousA }; } else if (fb * previousFb <= 0) { // sign change detected at near upper bound return new double[] { previousB, b }; } } } // no bracketing found throw new NoBracketingException(a, b, fa, fb); }
Compute the midpoint of two values.
Params:
  • a – first value.
  • b – second value.
Returns:the midpoint.
/** * Compute the midpoint of two values. * * @param a first value. * @param b second value. * @return the midpoint. */
public static double midpoint(double a, double b) { return (a + b) * 0.5; }
Check whether the interval bounds bracket a root. That is, if the values at the endpoints are not equal to zero, then the function takes opposite signs at the endpoints.
Params:
  • function – Function.
  • lower – Lower endpoint.
  • upper – Upper endpoint.
Throws:
Returns:true if the function values have opposite signs at the given points.
/** * Check whether the interval bounds bracket a root. That is, if the * values at the endpoints are not equal to zero, then the function takes * opposite signs at the endpoints. * * @param function Function. * @param lower Lower endpoint. * @param upper Upper endpoint. * @return {@code true} if the function values have opposite signs at the * given points. * @throws NullArgumentException if {@code function} is {@code null}. */
public static boolean isBracketing(UnivariateFunction function, final double lower, final double upper) throws NullArgumentException { if (function == null) { throw new NullArgumentException(LocalizedFormats.FUNCTION); } final double fLo = function.value(lower); final double fHi = function.value(upper); return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0); }
Check whether the arguments form a (strictly) increasing sequence.
Params:
  • start – First number.
  • mid – Second number.
  • end – Third number.
Returns:true if the arguments form an increasing sequence.
/** * Check whether the arguments form a (strictly) increasing sequence. * * @param start First number. * @param mid Second number. * @param end Third number. * @return {@code true} if the arguments form an increasing sequence. */
public static boolean isSequence(final double start, final double mid, final double end) { return (start < mid) && (mid < end); }
Check that the endpoints specify an interval.
Params:
  • lower – Lower endpoint.
  • upper – Upper endpoint.
Throws:
/** * Check that the endpoints specify an interval. * * @param lower Lower endpoint. * @param upper Upper endpoint. * @throws NumberIsTooLargeException if {@code lower >= upper}. */
public static void verifyInterval(final double lower, final double upper) throws NumberIsTooLargeException { if (lower >= upper) { throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, lower, upper, false); } }
Check that lower < initial < upper.
Params:
  • lower – Lower endpoint.
  • initial – Initial value.
  • upper – Upper endpoint.
Throws:
/** * Check that {@code lower < initial < upper}. * * @param lower Lower endpoint. * @param initial Initial value. * @param upper Upper endpoint. * @throws NumberIsTooLargeException if {@code lower >= initial} or * {@code initial >= upper}. */
public static void verifySequence(final double lower, final double initial, final double upper) throws NumberIsTooLargeException { verifyInterval(lower, initial); verifyInterval(initial, upper); }
Check that the endpoints specify an interval and the end points bracket a root.
Params:
  • function – Function.
  • lower – Lower endpoint.
  • upper – Upper endpoint.
Throws:
/** * Check that the endpoints specify an interval and the end points * bracket a root. * * @param function Function. * @param lower Lower endpoint. * @param upper Upper endpoint. * @throws NoBracketingException if the function has the same sign at the * endpoints. * @throws NullArgumentException if {@code function} is {@code null}. */
public static void verifyBracketing(UnivariateFunction function, final double lower, final double upper) throws NullArgumentException, NoBracketingException { if (function == null) { throw new NullArgumentException(LocalizedFormats.FUNCTION); } verifyInterval(lower, upper); if (!isBracketing(function, lower, upper)) { throw new NoBracketingException(lower, upper, function.value(lower), function.value(upper)); } } }