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