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