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    }