    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    import org.apache.commons.math.ConvergenceException;
    020    import org.apache.commons.math.FunctionEvaluationException;
    021    import org.apache.commons.math.MathRuntimeException;
    022    import org.apache.commons.math.analysis.UnivariateRealFunction;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.exception.NullArgumentException;
    025    import org.apache.commons.math.util.FastMath;
    026
    027    /**
    028     * Utility routines for {@link UnivariateRealSolver} objects.
    029     *
    030     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    031     */
    032    public class UnivariateRealSolverUtils {
    033
    034        /**
    035         * Default constructor.
    036         */
    037        private UnivariateRealSolverUtils() {
    038            super();
    039        }
    040
    041        /**
    042         * Convenience method to find a zero of a univariate real function.  A default
    043         * solver is used.
    044         *
    045         * @param f the function.
    046         * @param x0 the lower bound for the interval.
    047         * @param x1 the upper bound for the interval.
    048         * @return a value where the function is zero.
    049         * @throws ConvergenceException if the iteration count was exceeded
    050         * @throws FunctionEvaluationException if an error occurs evaluating the function
    051         * @throws IllegalArgumentException if f is null or the endpoints do not
    052         * specify a valid interval
    053         */
    054        public static double solve(UnivariateRealFunction f, double x0, double x1)
    055        throws ConvergenceException, FunctionEvaluationException {
    056            setup(f);
    057            return LazyHolder.FACTORY.newDefaultSolver().solve(f, 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 f the function
    065         * @param x0 the lower bound for the interval
    066         * @param x1 the upper bound for the interval
    067         * @param absoluteAccuracy the accuracy to be used by the solver
    068         * @return a value where the function is zero
    069         * @throws ConvergenceException if the iteration count is exceeded
    070         * @throws FunctionEvaluationException if an error occurs evaluating the function
    071         * @throws IllegalArgumentException if f is null, the endpoints do not
    072         * specify a valid interval, or the absoluteAccuracy is not valid for the
    073         * default solver
    074         */
    075        public static double solve(UnivariateRealFunction f, double x0, double x1,
    076                double absoluteAccuracy) throws ConvergenceException,
    077                FunctionEvaluationException {
    078
    079            setup(f);
    080            UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
    081            solver.setAbsoluteAccuracy(absoluteAccuracy);
    082            return solver.solve(f, x0, x1);
    083        }
    084
    085        /**
    086         * This method attempts to find two values a and b satisfying <ul>
    087        * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
    088         * <li> <code> f(a) * f(b) < 0 </code></li>
    089         * </ul>
    090         * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
    091         * and <code>b</code> bracket a root of f.
    092         * <p>
    093         * The algorithm starts by setting
    094         * <code>a := initial -1; b := initial +1,</code> examines the value of the
    095         * function at <code>a</code> and <code>b</code> and keeps moving
    096         * the endpoints out by one unit each time through a loop that terminates
    097         * when one of the following happens: <ul>
    098         * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
    099         * <li> <code> a = lower </code> and <code> b = upper</code>
    100         * -- ConvergenceException </li>
    101         * <li> <code> Integer.MAX_VALUE</code> iterations elapse
    102         * -- ConvergenceException </li>
    103         * </ul></p>
    104         * <p>
    105         * <strong>Note: </strong> this method can take
    106         * <code>Integer.MAX_VALUE</code> iterations to throw a
    107         * <code>ConvergenceException.</code>  Unless you are confident that there
    108         * is a root between <code>lowerBound</code> and <code>upperBound</code>
    109         * near <code>initial,</code> it is better to use
    110         * {@link #bracket(UnivariateRealFunction, double, double, double, int)},
    111         * explicitly specifying the maximum number of iterations.</p>
    112         *
    113         * @param function the function
    114         * @param initial initial midpoint of interval being expanded to
    115         * bracket a root
    116         * @param lowerBound lower bound (a is never lower than this value)
    117         * @param upperBound upper bound (b never is greater than this
    118         * value)
    119         * @return a two element array holding {a, b}
    120         * @throws ConvergenceException if a root can not be bracketted
    121         * @throws FunctionEvaluationException if an error occurs evaluating the function
    122         * @throws IllegalArgumentException if function is null, maximumIterations
    123         * is not positive, or initial is not between lowerBound and upperBound
    124         */
    125        public static double[] bracket(UnivariateRealFunction function,
    126                double initial, double lowerBound, double upperBound)
    127        throws ConvergenceException, FunctionEvaluationException {
    128            return bracket( function, initial, lowerBound, upperBound,
    129                Integer.MAX_VALUE ) ;
    130        }
    131
    132         /**
    133         * This method attempts to find two values a and b satisfying <ul>
    134         * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
    135         * <li> <code> f(a) * f(b) <= 0 </code> </li>
    136         * </ul>
    137         * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
    138         * and <code>b</code> bracket a root of f.
    139         * <p>
    140         * The algorithm starts by setting
    141         * <code>a := initial -1; b := initial +1,</code> examines the value of the
    142         * function at <code>a</code> and <code>b</code> and keeps moving
    143         * the endpoints out by one unit each time through a loop that terminates
    144         * when one of the following happens: <ul>
    145         * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
    146         * <li> <code> a = lower </code> and <code> b = upper</code>
    147         * -- ConvergenceException </li>
    148         * <li> <code> maximumIterations</code> iterations elapse
    149         * -- ConvergenceException </li></ul></p>
    150         *
    151         * @param function the function
    152         * @param initial initial midpoint of interval being expanded to
    153         * bracket a root
    154         * @param lowerBound lower bound (a is never lower than this value)
    155         * @param upperBound upper bound (b never is greater than this
    156         * value)
    157         * @param maximumIterations maximum number of iterations to perform
    158         * @return a two element array holding {a, b}.
    159         * @throws ConvergenceException if the algorithm fails to find a and b
    160         * satisfying the desired conditions
    161         * @throws FunctionEvaluationException if an error occurs evaluating the function
    162         * @throws IllegalArgumentException if function is null, maximumIterations
    163         * is not positive, or initial is not between lowerBound and upperBound
    164         */
    165        public static double[] bracket(UnivariateRealFunction function,
    166                double initial, double lowerBound, double upperBound,
    167                int maximumIterations) throws ConvergenceException,
    168                FunctionEvaluationException {
    169
    170            if (function == null) {
    171                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    172            }
    173            if (maximumIterations <= 0)  {
    174                throw MathRuntimeException.createIllegalArgumentException(
    175                      LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    176            }
    177            if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
    178                throw MathRuntimeException.createIllegalArgumentException(
    179                      LocalizedFormats.INVALID_BRACKETING_PARAMETERS,
    180                      lowerBound, initial, upperBound);
    181            }
    182            double a = initial;
    183            double b = initial;
    184            double fa;
    185            double fb;
    186            int numIterations = 0 ;
    187
    188            do {
    189                a = FastMath.max(a - 1.0, lowerBound);
    190                b = FastMath.min(b + 1.0, upperBound);
    191                fa = function.value(a);
    192
    193                fb = function.value(b);
    194                numIterations++ ;
    195            } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
    196                    ((a > lowerBound) || (b < upperBound)));
    197
    198            if (fa * fb > 0.0 ) {
    199                throw new ConvergenceException(
    200                          LocalizedFormats.FAILED_BRACKETING,
    201                          numIterations, maximumIterations, initial,
    202                          lowerBound, upperBound, a, b, fa, fb);
    203            }
    204
    205            return new double[]{a, b};
    206        }
    207
    208        /**
    209         * Compute the midpoint of two values.
    210         *
    211         * @param a first value.
    212         * @param b second value.
    213         * @return the midpoint.
    214         */
    215        public static double midpoint(double a, double b) {
    216            return (a + b) * .5;
    217        }
    218
    219        /**
    220         * Checks to see if f is null, throwing IllegalArgumentException if so.
    221         * @param f  input function
    222         * @throws IllegalArgumentException if f is null
    223         */
    224        private static void setup(UnivariateRealFunction f) {
    225            if (f == null) {
    226                throw new NullArgumentException(LocalizedFormats.FUNCTION);
    227            }
    228        }
    229
    230        // CHECKSTYLE: stop HideUtilityClassConstructor
    231        /** Holder for the factory.
    232         * <p>We use here the Initialization On Demand Holder Idiom.</p>
    233         */
    234        private static class LazyHolder {
    235            /** Cached solver factory */
    236            private static final UnivariateRealSolverFactory FACTORY = UnivariateRealSolverFactory.newInstance();
    237        }
    238        // CHECKSTYLE: resume HideUtilityClassConstructor
    239
    240    }