    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.analysis.solvers;
    018
    019
    020    import org.apache.commons.math.FunctionEvaluationException;
    021    import org.apache.commons.math.MathRuntimeException;
    022    import org.apache.commons.math.MaxIterationsExceededException;
    023    import org.apache.commons.math.analysis.UnivariateRealFunction;
    024    import org.apache.commons.math.exception.util.LocalizedFormats;
    025    import org.apache.commons.math.util.FastMath;
    026
    027    /**
    028     * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
    029     * Brent algorithm</a> for  finding zeros of real univariate functions.
    030     * <p>
    031     * The function should be continuous but not necessarily smooth.</p>
    032     *
    033     * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
    034     */
    035    public class BrentSolver extends UnivariateRealSolverImpl {
    036
    037        /**
    038         * Default absolute accuracy
    039         * @since 2.1
    040         */
    041        public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
    042
    043        /** Default maximum number of iterations
    044         * @since 2.1
    045         */
    046        public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
    047
    048        /** Serializable version identifier */
    049        private static final long serialVersionUID = 7694577816772532779L;
    050
    051        /**
    052         * Construct a solver for the given function.
    053         *
    054         * @param f function to solve.
    055         * @deprecated as of 2.0 the function to solve is passed as an argument
    056         * to the {@link #solve(UnivariateRealFunction, double, double)} or
    057         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
    058         * method.
    059         */
    060        @Deprecated
    061        public BrentSolver(UnivariateRealFunction f) {
    062            super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
    063        }
    064
    065        /**
    066         * Construct a solver with default properties.
    067         * @deprecated in 2.2 (to be removed in 3.0).
    068         */
    069        @Deprecated
    070        public BrentSolver() {
    071            super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
    072        }
    073
    074        /**
    075         * Construct a solver with the given absolute accuracy.
    076         *
    077         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
    078         * @since 2.1
    079         */
    080        public BrentSolver(double absoluteAccuracy) {
    081            super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
    082        }
    083
    084        /**
    085         * Contstruct a solver with the given maximum iterations and absolute accuracy.
    086         *
    087         * @param maximumIterations maximum number of iterations
    088         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
    089         * @since 2.1
    090         */
    091        public BrentSolver(int maximumIterations, double absoluteAccuracy) {
    092            super(maximumIterations, absoluteAccuracy);
    093        }
    094
    095        /** {@inheritDoc} */
    096        @Deprecated
    097        public double solve(double min, double max)
    098            throws MaxIterationsExceededException, FunctionEvaluationException {
    099            return solve(f, min, max);
    100        }
    101
    102        /** {@inheritDoc} */
    103        @Deprecated
    104        public double solve(double min, double max, double initial)
    105            throws MaxIterationsExceededException, FunctionEvaluationException {
    106            return solve(f, min, max, initial);
    107        }
    108
    109        /**
    110         * Find a zero in the given interval with an initial guess.
    111         * <p>Throws <code>IllegalArgumentException</code> if the values of the
    112         * function at the three points have the same sign (note that it is
    113         * allowed to have endpoints with the same sign if the initial point has
    114         * opposite sign function-wise).</p>
    115         *
    116         * @param f function to solve.
    117         * @param min the lower bound for the interval.
    118         * @param max the upper bound for the interval.
    119         * @param initial the start value to use (must be set to min if no
    120         * initial point is known).
    121         * @return the value where the function is zero
    122         * @throws MaxIterationsExceededException the maximum iteration count is exceeded
    123         * @throws FunctionEvaluationException if an error occurs evaluating  the function
    124         * @throws IllegalArgumentException if initial is not between min and max
    125         * (even if it <em>is</em> a root)
    126         * @deprecated in 2.2 (to be removed in 3.0).
    127         */
    128        @Deprecated
    129        public double solve(final UnivariateRealFunction f,
    130                            final double min, final double max, final double initial)
    131            throws MaxIterationsExceededException, FunctionEvaluationException {
    132
    133            clearResult();
    134            if ((initial < min) || (initial > max)) {
    135                throw MathRuntimeException.createIllegalArgumentException(
    136                      LocalizedFormats.INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS,
    137                      min, initial, max);
    138            }
    139
    140            // return the initial guess if it is good enough
    141            double yInitial = f.value(initial);
    142            if (FastMath.abs(yInitial) <= functionValueAccuracy) {
    143                setResult(initial, 0);
    144                return result;
    145            }
    146
    147            // return the first endpoint if it is good enough
    148            double yMin = f.value(min);
    149            if (FastMath.abs(yMin) <= functionValueAccuracy) {
    150                setResult(min, 0);
    151                return result;
    152            }
    153
    154            // reduce interval if min and initial bracket the root
    155            if (yInitial * yMin < 0) {
    156                return solve(f, min, yMin, initial, yInitial, min, yMin);
    157            }
    158
    159            // return the second endpoint if it is good enough
    160            double yMax = f.value(max);
    161            if (FastMath.abs(yMax) <= functionValueAccuracy) {
    162                setResult(max, 0);
    163                return result;
    164            }
    165
    166            // reduce interval if initial and max bracket the root
    167            if (yInitial * yMax < 0) {
    168                return solve(f, initial, yInitial, max, yMax, initial, yInitial);
    169            }
    170
    171            throw MathRuntimeException.createIllegalArgumentException(
    172                  LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
    173
    174        }
    175
    176        /**
    177         * Find a zero in the given interval with an initial guess.
    178         * <p>Throws <code>IllegalArgumentException</code> if the values of the
    179         * function at the three points have the same sign (note that it is
    180         * allowed to have endpoints with the same sign if the initial point has
    181         * opposite sign function-wise).</p>
    182         *
    183         * @param f function to solve.
    184         * @param min the lower bound for the interval.
    185         * @param max the upper bound for the interval.
    186         * @param initial the start value to use (must be set to min if no
    187         * initial point is known).
    188         * @param maxEval Maximum number of evaluations.
    189         * @return the value where the function is zero
    190         * @throws MaxIterationsExceededException the maximum iteration count is exceeded
    191         * @throws FunctionEvaluationException if an error occurs evaluating  the function
    192         * @throws IllegalArgumentException if initial is not between min and max
    193         * (even if it <em>is</em> a root)
    194         */
    195        @Override
    196        public double solve(int maxEval, final UnivariateRealFunction f,
    197                            final double min, final double max, final double initial)
    198            throws MaxIterationsExceededException, FunctionEvaluationException {
    199            setMaximalIterationCount(maxEval);
    200            return solve(f, min, max, initial);
    201        }
    202
    203        /**
    204         * Find a zero in the given interval.
    205         * <p>
    206         * Requires that the values of the function at the endpoints have opposite
    207         * signs. An <code>IllegalArgumentException</code> is thrown if this is not
    208         * the case.</p>
    209         *
    210         * @param f the function to solve
    211         * @param min the lower bound for the interval.
    212         * @param max the upper bound for the interval.
    213         * @return the value where the function is zero
    214         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    215         * @throws FunctionEvaluationException if an error occurs evaluating the function
    216         * @throws IllegalArgumentException if min is not less than max or the
    217         * signs of the values of the function at the endpoints are not opposites
    218         * @deprecated in 2.2 (to be removed in 3.0).
    219         */
    220        @Deprecated
    221        public double solve(final UnivariateRealFunction f,
    222                            final double min, final double max)
    223            throws MaxIterationsExceededException, FunctionEvaluationException {
    224
    225            clearResult();
    226            verifyInterval(min, max);
    227
    228            double ret = Double.NaN;
    229
    230            double yMin = f.value(min);
    231            double yMax = f.value(max);
    232
    233            // Verify bracketing
    234            double sign = yMin * yMax;
    235            if (sign > 0) {
    236                // check if either value is close to a zero
    237                if (FastMath.abs(yMin) <= functionValueAccuracy) {
    238                    setResult(min, 0);
    239                    ret = min;
    240                } else if (FastMath.abs(yMax) <= functionValueAccuracy) {
    241                    setResult(max, 0);
    242                    ret = max;
    243                } else {
    244                    // neither value is close to zero and min and max do not bracket root.
    245                    throw MathRuntimeException.createIllegalArgumentException(
    246                            LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
    247                }
    248            } else if (sign < 0){
    249                // solve using only the first endpoint as initial guess
    250                ret = solve(f, min, yMin, max, yMax, min, yMin);
    251            } else {
    252                // either min or max is a root
    253                if (yMin == 0.0) {
    254                    ret = min;
    255                } else {
    256                    ret = max;
    257                }
    258            }
    259
    260            return ret;
    261        }
    262
    263        /**
    264         * Find a zero in the given interval.
    265         * <p>
    266         * Requires that the values of the function at the endpoints have opposite
    267         * signs. An <code>IllegalArgumentException</code> is thrown if this is not
    268         * the case.</p>
    269         *
    270         * @param f the function to solve
    271         * @param min the lower bound for the interval.
    272         * @param max the upper bound for the interval.
    273         * @param maxEval Maximum number of evaluations.
    274         * @return the value where the function is zero
    275         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    276         * @throws FunctionEvaluationException if an error occurs evaluating the function
    277         * @throws IllegalArgumentException if min is not less than max or the
    278         * signs of the values of the function at the endpoints are not opposites
    279         */
    280        @Override
    281        public double solve(int maxEval, final UnivariateRealFunction f,
    282                            final double min, final double max)
    283            throws MaxIterationsExceededException, FunctionEvaluationException {
    284            setMaximalIterationCount(maxEval);
    285            return solve(f, min, max);
    286        }
    287
    288        /**
    289         * Find a zero starting search according to the three provided points.
    290         * @param f the function to solve
    291         * @param x0 old approximation for the root
    292         * @param y0 function value at the approximation for the root
    293         * @param x1 last calculated approximation for the root
    294         * @param y1 function value at the last calculated approximation
    295         * for the root
    296         * @param x2 bracket point (must be set to x0 if no bracket point is
    297         * known, this will force starting with linear interpolation)
    298         * @param y2 function value at the bracket point.
    299         * @return the value where the function is zero
    300         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    301         * @throws FunctionEvaluationException if an error occurs evaluating the function
    302         */
    303        private double solve(final UnivariateRealFunction f,
    304                             double x0, double y0,
    305                             double x1, double y1,
    306                             double x2, double y2)
    307        throws MaxIterationsExceededException, FunctionEvaluationException {
    308
    309            double delta = x1 - x0;
    310            double oldDelta = delta;
    311
    312            int i = 0;
    313            while (i < maximalIterationCount) {
    314                if (FastMath.abs(y2) < FastMath.abs(y1)) {
    315                    // use the bracket point if is better than last approximation
    316                    x0 = x1;
    317                    x1 = x2;
    318                    x2 = x0;
    319                    y0 = y1;
    320                    y1 = y2;
    321                    y2 = y0;
    322                }
    323                if (FastMath.abs(y1) <= functionValueAccuracy) {
    324                    // Avoid division by very small values. Assume
    325                    // the iteration has converged (the problem may
    326                    // still be ill conditioned)
    327                    setResult(x1, i);
    328                    return result;
    329                }
    330                double dx = x2 - x1;
    331                double tolerance =
    332                    FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy);
    333                if (FastMath.abs(dx) <= tolerance) {
    334                    setResult(x1, i);
    335                    return result;
    336                }
    337                if ((FastMath.abs(oldDelta) < tolerance) ||
    338                        (FastMath.abs(y0) <= FastMath.abs(y1))) {
    339                    // Force bisection.
    340                    delta = 0.5 * dx;
    341                    oldDelta = delta;
    342                } else {
    343                    double r3 = y1 / y0;
    344                    double p;
    345                    double p1;
    346                    // the equality test (x0 == x2) is intentional,
    347                    // it is part of the original Brent's method,
    348                    // it should NOT be replaced by proximity test
    349                    if (x0 == x2) {
    350                        // Linear interpolation.
    351                        p = dx * r3;
    352                        p1 = 1.0 - r3;
    353                    } else {
    354                        // Inverse quadratic interpolation.
    355                        double r1 = y0 / y2;
    356                        double r2 = y1 / y2;
    357                        p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
    358                        p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
    359                    }
    360                    if (p > 0.0) {
    361                        p1 = -p1;
    362                    } else {
    363                        p = -p;
    364                    }
    365                    if (2.0 * p >= 1.5 * dx * p1 - FastMath.abs(tolerance * p1) ||
    366                            p >= FastMath.abs(0.5 * oldDelta * p1)) {
    367                        // Inverse quadratic interpolation gives a value
    368                        // in the wrong direction, or progress is slow.
    369                        // Fall back to bisection.
    370                        delta = 0.5 * dx;
    371                        oldDelta = delta;
    372                    } else {
    373                        oldDelta = delta;
    374                        delta = p / p1;
    375                    }
    376                }
    377                // Save old X1, Y1
    378                x0 = x1;
    379                y0 = y1;
    380                // Compute new X1, Y1
    381                if (FastMath.abs(delta) > tolerance) {
    382                    x1 = x1 + delta;
    383                } else if (dx > 0.0) {
    384                    x1 = x1 + 0.5 * tolerance;
    385                } else if (dx <= 0.0) {
    386                    x1 = x1 - 0.5 * tolerance;
    387                }
    388                y1 = f.value(x1);
    389                if ((y1 > 0) == (y2 > 0)) {
    390                    x2 = x0;
    391                    y2 = y0;
    392                    delta = x1 - x0;
    393                    oldDelta = delta;
    394                }
    395                i++;
    396            }
    397            throw new MaxIterationsExceededException(maximalIterationCount);
    398        }
    399    }