View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.analysis.solvers;
18  
19  import org.apache.commons.math.FunctionEvaluationException;
20  import org.apache.commons.math.ConvergenceException;
21  import org.apache.commons.math.MathRuntimeException;
22  import org.apache.commons.math.analysis.UnivariateRealFunction;
23  
24  /**
25   * Utility routines for {@link UnivariateRealSolver} objects.
26   * 
27   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
28   */
29  public class UnivariateRealSolverUtils {
30      /**
31       * Default constructor.
32       */
33      private UnivariateRealSolverUtils() {
34          super();
35      }
36      
37      /**
38       * Convenience method to find a zero of a univariate real function.  A default
39       * solver is used. 
40       * 
41       * @param f the function.
42       * @param x0 the lower bound for the interval.
43       * @param x1 the upper bound for the interval.
44       * @return a value where the function is zero.
45       * @throws ConvergenceException if the iteration count was exceeded
46       * @throws FunctionEvaluationException if an error occurs evaluating
47       * the function
48       * @throws IllegalArgumentException if f is null or the endpoints do not
49       * specify a valid interval
50       */
51      public static double solve(UnivariateRealFunction f, double x0, double x1)
52      throws ConvergenceException, FunctionEvaluationException {
53          setup(f);
54          return LazyHolder.FACTORY.newDefaultSolver().solve(f, x0, x1);
55      }
56  
57      /**
58       * Convenience method to find a zero of a univariate real function.  A default
59       * solver is used. 
60       * 
61       * @param f the function
62       * @param x0 the lower bound for the interval
63       * @param x1 the upper bound for the interval
64       * @param absoluteAccuracy the accuracy to be used by the solver
65       * @return a value where the function is zero
66       * @throws ConvergenceException if the iteration count is exceeded
67       * @throws FunctionEvaluationException if an error occurs evaluating the
68       * function
69       * @throws IllegalArgumentException if f is null, the endpoints do not 
70       * specify a valid interval, or the absoluteAccuracy is not valid for the
71       * default solver
72       */
73      public static double solve(UnivariateRealFunction f, double x0, double x1,
74              double absoluteAccuracy) throws ConvergenceException, 
75              FunctionEvaluationException {    
76         
77          setup(f);
78          UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
79          solver.setAbsoluteAccuracy(absoluteAccuracy);
80          return solver.solve(f, x0, x1);
81      }
82  
83      /**
84       * This method attempts to find two values a and b satisfying <ul>
85      * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
86       * <li> <code> f(a) * f(b) < 0 </code></li>
87       * </ul>
88       * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
89       * and <code>b</code> bracket a root of f.
90       * <p>
91       * The algorithm starts by setting 
92       * <code>a := initial -1; b := initial +1,</code> examines the value of the
93       * function at <code>a</code> and <code>b</code> and keeps moving
94       * the endpoints out by one unit each time through a loop that terminates 
95       * when one of the following happens: <ul>
96       * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
97       * <li> <code> a = lower </code> and <code> b = upper</code> 
98       * -- ConvergenceException </li>
99       * <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 }