    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.math3.analysis.solvers;
    018
    019    import org.apache.commons.math3.analysis.UnivariateFunction;
    020    import org.apache.commons.math3.exception.NoBracketingException;
    021    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    022    import org.apache.commons.math3.exception.NullArgumentException;
    023    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.util.FastMath;
    026
    027    /**
    028     * Utility routines for {@link UnivariateSolver} objects.
    029     *
    030     * @version $Id: UnivariateSolverUtils.java 1400850 2012-10-22 11:57:17Z erans $
    031     */
    032    public class UnivariateSolverUtils {
    033        /**
    034         * Class contains only static methods.
    035         */
    036        private UnivariateSolverUtils() {}
    037
    038        /**
    039         * Convenience method to find a zero of a univariate real function.  A default
    040         * solver is used.
    041         *
    042         * @param function Function.
    043         * @param x0 Lower bound for the interval.
    044         * @param x1 Upper bound for the interval.
    045         * @return a value where the function is zero.
    046         * @throws NoBracketingException if the function has the same sign at the
    047         * endpoints.
    048         * @throws NullArgumentException if {@code function} is {@code null}.
    049         */
    050        public static double solve(UnivariateFunction function, double x0, double x1)
    051            throws NullArgumentException,
    052                   NoBracketingException {
    053            if (function == null) {
    054                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    055            }
    056            final UnivariateSolver solver = new BrentSolver();
    057            return solver.solve(Integer.MAX_VALUE, function, x0, x1);
    058        }
    059
    060        /**
    061         * Convenience method to find a zero of a univariate real function.  A default
    062         * solver is used.
    063         *
    064         * @param function Function.
    065         * @param x0 Lower bound for the interval.
    066         * @param x1 Upper bound for the interval.
    067         * @param absoluteAccuracy Accuracy to be used by the solver.
    068         * @return a value where the function is zero.
    069         * @throws NoBracketingException if the function has the same sign at the
    070         * endpoints.
    071         * @throws NullArgumentException if {@code function} is {@code null}.
    072         */
    073        public static double solve(UnivariateFunction function,
    074                                   double x0, double x1,
    075                                   double absoluteAccuracy)
    076            throws NullArgumentException,
    077                   NoBracketingException {
    078            if (function == null) {
    079                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    080            }
    081            final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
    082            return solver.solve(Integer.MAX_VALUE, function, x0, x1);
    083        }
    084
    085        /** Force a root found by a non-bracketing solver to lie on a specified side,
    086         * as if the solver was a bracketing one.
    087         * @param maxEval maximal number of new evaluations of the function
    088         * (evaluations already done for finding the root should have already been subtracted
    089         * from this number)
    090         * @param f function to solve
    091         * @param bracketing bracketing solver to use for shifting the root
    092         * @param baseRoot original root found by a previous non-bracketing solver
    093         * @param min minimal bound of the search interval
    094         * @param max maximal bound of the search interval
    095         * @param allowedSolution the kind of solutions that the root-finding algorithm may
    096         * accept as solutions.
    097         * @return a root approximation, on the specified side of the exact root
    098         * @throws NoBracketingException if the function has the same sign at the
    099         * endpoints.
    100         */
    101        public static double forceSide(final int maxEval, final UnivariateFunction f,
    102                                       final BracketedUnivariateSolver<UnivariateFunction> bracketing,
    103                                       final double baseRoot, final double min, final double max,
    104                                       final AllowedSolution allowedSolution)
    105            throws NoBracketingException {
    106
    107            if (allowedSolution == AllowedSolution.ANY_SIDE) {
    108                // no further bracketing required
    109                return baseRoot;
    110            }
    111
    112            // find a very small interval bracketing the root
    113            final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
    114                                             FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
    115            double xLo        = FastMath.max(min, baseRoot - step);
    116            double fLo        = f.value(xLo);
    117            double xHi        = FastMath.min(max, baseRoot + step);
    118            double fHi        = f.value(xHi);
    119            int remainingEval = maxEval - 2;
    120            while (remainingEval > 0) {
    121
    122                if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
    123                    // compute the root on the selected side
    124                    return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
    125                }
    126
    127                // try increasing the interval
    128                boolean changeLo = false;
    129                boolean changeHi = false;
    130                if (fLo < fHi) {
    131                    // increasing function
    132                    if (fLo >= 0) {
    133                        changeLo = true;
    134                    } else {
    135                        changeHi = true;
    136                    }
    137                } else if (fLo > fHi) {
    138                    // decreasing function
    139                    if (fLo <= 0) {
    140                        changeLo = true;
    141                    } else {
    142                        changeHi = true;
    143                    }
    144                } else {
    145                    // unknown variation
    146                    changeLo = true;
    147                    changeHi = true;
    148                }
    149
    150                // update the lower bound
    151                if (changeLo) {
    152                    xLo = FastMath.max(min, xLo - step);
    153                    fLo  = f.value(xLo);
    154                    remainingEval--;
    155                }
    156
    157                // update the higher bound
    158                if (changeHi) {
    159                    xHi = FastMath.min(max, xHi + step);
    160                    fHi  = f.value(xHi);
    161                    remainingEval--;
    162                }
    163
    164            }
    165
    166            throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
    167                                            xLo, xHi, fLo, fHi,
    168                                            maxEval - remainingEval, maxEval, baseRoot,
    169                                            min, max);
    170
    171        }
    172
    173        /**
    174         * This method attempts to find two values a and b satisfying <ul>
    175         * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
    176         * <li> <code> f(a) * f(b) < 0 </code></li>
    177         * </ul>
    178         * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
    179         * and <code>b</code> bracket a root of f.
    180         * <p>
    181         * The algorithm starts by setting
    182         * <code>a := initial -1; b := initial +1,</code> examines the value of the
    183         * function at <code>a</code> and <code>b</code> and keeps moving
    184         * the endpoints out by one unit each time through a loop that terminates
    185         * when one of the following happens: <ul>
    186         * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
    187         * <li> <code> a = lower </code> and <code> b = upper</code>
    188         * -- NoBracketingException </li>
    189         * <li> <code> Integer.MAX_VALUE</code> iterations elapse
    190         * -- NoBracketingException </li>
    191         * </ul></p>
    192         * <p>
    193         * <strong>Note: </strong> this method can take
    194         * <code>Integer.MAX_VALUE</code> iterations to throw a
    195         * <code>ConvergenceException.</code>  Unless you are confident that there
    196         * is a root between <code>lowerBound</code> and <code>upperBound</code>
    197         * near <code>initial,</code> it is better to use
    198         * {@link #bracket(UnivariateFunction, double, double, double, int)},
    199         * explicitly specifying the maximum number of iterations.</p>
    200         *
    201         * @param function Function.
    202         * @param initial Initial midpoint of interval being expanded to
    203         * bracket a root.
    204         * @param lowerBound Lower bound (a is never lower than this value)
    205         * @param upperBound Upper bound (b never is greater than this
    206         * value).
    207         * @return a two-element array holding a and b.
    208         * @throws NoBracketingException if a root cannot be bracketted.
    209         * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
    210         * @throws NullArgumentException if {@code function} is {@code null}.
    211         */
    212        public static double[] bracket(UnivariateFunction function,
    213                                       double initial,
    214                                       double lowerBound, double upperBound)
    215            throws NullArgumentException,
    216                   NotStrictlyPositiveException,
    217                   NoBracketingException {
    218            return bracket(function, initial, lowerBound, upperBound, Integer.MAX_VALUE);
    219        }
    220
    221         /**
    222         * This method attempts to find two values a and b satisfying <ul>
    223         * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
    224         * <li> <code> f(a) * f(b) <= 0 </code> </li>
    225         * </ul>
    226         * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
    227         * and <code>b</code> bracket a root of f.
    228         * <p>
    229         * The algorithm starts by setting
    230         * <code>a := initial -1; b := initial +1,</code> examines the value of the
    231         * function at <code>a</code> and <code>b</code> and keeps moving
    232         * the endpoints out by one unit each time through a loop that terminates
    233         * when one of the following happens: <ul>
    234         * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
    235         * <li> <code> a = lower </code> and <code> b = upper</code>
    236         * -- NoBracketingException </li>
    237         * <li> <code> maximumIterations</code> iterations elapse
    238         * -- NoBracketingException </li></ul></p>
    239         *
    240         * @param function Function.
    241         * @param initial Initial midpoint of interval being expanded to
    242         * bracket a root.
    243         * @param lowerBound Lower bound (a is never lower than this value).
    244         * @param upperBound Upper bound (b never is greater than this
    245         * value).
    246         * @param maximumIterations Maximum number of iterations to perform
    247         * @return a two element array holding a and b.
    248         * @throws NoBracketingException if the algorithm fails to find a and b
    249         * satisfying the desired conditions.
    250         * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
    251         * @throws NullArgumentException if {@code function} is {@code null}.
    252         */
    253        public static double[] bracket(UnivariateFunction function,
    254                                       double initial,
    255                                       double lowerBound, double upperBound,
    256                                       int maximumIterations)
    257            throws NullArgumentException,
    258                   NotStrictlyPositiveException,
    259                   NoBracketingException {
    260            if (function == null) {
    261                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    262            }
    263            if (maximumIterations <= 0)  {
    264                throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    265            }
    266            verifySequence(lowerBound, initial, upperBound);
    267
    268            double a = initial;
    269            double b = initial;
    270            double fa;
    271            double fb;
    272            int numIterations = 0;
    273
    274            do {
    275                a = FastMath.max(a - 1.0, lowerBound);
    276                b = FastMath.min(b + 1.0, upperBound);
    277                fa = function.value(a);
    278
    279                fb = function.value(b);
    280                ++numIterations;
    281            } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
    282                    ((a > lowerBound) || (b < upperBound)));
    283
    284            if (fa * fb > 0.0) {
    285                throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
    286                                                a, b, fa, fb,
    287                                                numIterations, maximumIterations, initial,
    288                                                lowerBound, upperBound);
    289            }
    290
    291            return new double[] {a, b};
    292        }
    293
    294        /**
    295         * Compute the midpoint of two values.
    296         *
    297         * @param a first value.
    298         * @param b second value.
    299         * @return the midpoint.
    300         */
    301        public static double midpoint(double a, double b) {
    302            return (a + b) * 0.5;
    303        }
    304
    305        /**
    306         * Check whether the interval bounds bracket a root. That is, if the
    307         * values at the endpoints are not equal to zero, then the function takes
    308         * opposite signs at the endpoints.
    309         *
    310         * @param function Function.
    311         * @param lower Lower endpoint.
    312         * @param upper Upper endpoint.
    313         * @return {@code true} if the function values have opposite signs at the
    314         * given points.
    315         * @throws NullArgumentException if {@code function} is {@code null}.
    316         */
    317        public static boolean isBracketing(UnivariateFunction function,
    318                                           final double lower,
    319                                           final double upper)
    320            throws NullArgumentException {
    321            if (function == null) {
    322                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    323            }
    324            final double fLo = function.value(lower);
    325            final double fHi = function.value(upper);
    326            return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
    327        }
    328
    329        /**
    330         * Check whether the arguments form a (strictly) increasing sequence.
    331         *
    332         * @param start First number.
    333         * @param mid Second number.
    334         * @param end Third number.
    335         * @return {@code true} if the arguments form an increasing sequence.
    336         */
    337        public static boolean isSequence(final double start,
    338                                         final double mid,
    339                                         final double end) {
    340            return (start < mid) && (mid < end);
    341        }
    342
    343        /**
    344         * Check that the endpoints specify an interval.
    345         *
    346         * @param lower Lower endpoint.
    347         * @param upper Upper endpoint.
    348         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    349         */
    350        public static void verifyInterval(final double lower,
    351                                          final double upper)
    352            throws NumberIsTooLargeException {
    353            if (lower >= upper) {
    354                throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
    355                                                    lower, upper, false);
    356            }
    357        }
    358
    359        /**
    360         * Check that {@code lower < initial < upper}.
    361         *
    362         * @param lower Lower endpoint.
    363         * @param initial Initial value.
    364         * @param upper Upper endpoint.
    365         * @throws NumberIsTooLargeException if {@code lower >= initial} or
    366         * {@code initial >= upper}.
    367         */
    368        public static void verifySequence(final double lower,
    369                                          final double initial,
    370                                          final double upper)
    371            throws NumberIsTooLargeException {
    372            verifyInterval(lower, initial);
    373            verifyInterval(initial, upper);
    374        }
    375
    376        /**
    377         * Check that the endpoints specify an interval and the end points
    378         * bracket a root.
    379         *
    380         * @param function Function.
    381         * @param lower Lower endpoint.
    382         * @param upper Upper endpoint.
    383         * @throws NoBracketingException if the function has the same sign at the
    384         * endpoints.
    385         * @throws NullArgumentException if {@code function} is {@code null}.
    386         */
    387        public static void verifyBracketing(UnivariateFunction function,
    388                                            final double lower,
    389                                            final double upper)
    390            throws NullArgumentException,
    391                   NoBracketingException {
    392            if (function == null) {
    393                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    394            }
    395            verifyInterval(lower, upper);
    396            if (!isBracketing(function, lower, upper)) {
    397                throw new NoBracketingException(lower, upper,
    398                                                function.value(lower),
    399                                                function.value(upper));
    400            }
    401        }
    402    }