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.math3.analysis.solvers;
18  
19  import org.apache.commons.math3.analysis.UnivariateFunction;
20  import org.apache.commons.math3.exception.NoBracketingException;
21  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
22  import org.apache.commons.math3.exception.NullArgumentException;
23  import org.apache.commons.math3.exception.NumberIsTooLargeException;
24  import org.apache.commons.math3.exception.util.LocalizedFormats;
25  import org.apache.commons.math3.util.FastMath;
26  
27  /**
28   * Utility routines for {@link UnivariateSolver} objects.
29   *
30   * @version $Id: UnivariateSolverUtils.java 1400850 2012-10-22 11:57:17Z erans $
31   */
32  public class UnivariateSolverUtils {
33      /**
34       * Class contains only static methods.
35       */
36      private UnivariateSolverUtils() {}
37  
38      /**
39       * Convenience method to find a zero of a univariate real function.  A default
40       * solver is used.
41       *
42       * @param function Function.
43       * @param x0 Lower bound for the interval.
44       * @param x1 Upper bound for the interval.
45       * @return a value where the function is zero.
46       * @throws NoBracketingException if the function has the same sign at the
47       * endpoints.
48       * @throws NullArgumentException if {@code function} is {@code null}.
49       */
50      public static double solve(UnivariateFunction function, double x0, double x1)
51          throws NullArgumentException,
52                 NoBracketingException {
53          if (function == null) {
54              throw new NullArgumentException(LocalizedFormats.FUNCTION);
55          }
56          final UnivariateSolver solver = new BrentSolver();
57          return solver.solve(Integer.MAX_VALUE, function, x0, x1);
58      }
59  
60      /**
61       * Convenience method to find a zero of a univariate real function.  A default
62       * solver is used.
63       *
64       * @param function Function.
65       * @param x0 Lower bound for the interval.
66       * @param x1 Upper bound for the interval.
67       * @param absoluteAccuracy Accuracy to be used by the solver.
68       * @return a value where the function is zero.
69       * @throws NoBracketingException if the function has the same sign at the
70       * endpoints.
71       * @throws NullArgumentException if {@code function} is {@code null}.
72       */
73      public static double solve(UnivariateFunction function,
74                                 double x0, double x1,
75                                 double absoluteAccuracy)
76          throws NullArgumentException,
77                 NoBracketingException {
78          if (function == null) {
79              throw new NullArgumentException(LocalizedFormats.FUNCTION);
80          }
81          final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
82          return solver.solve(Integer.MAX_VALUE, function, x0, x1);
83      }
84  
85      /** Force a root found by a non-bracketing solver to lie on a specified side,
86       * as if the solver was a bracketing one.
87       * @param maxEval maximal number of new evaluations of the function
88       * (evaluations already done for finding the root should have already been subtracted
89       * from this number)
90       * @param f function to solve
91       * @param bracketing bracketing solver to use for shifting the root
92       * @param baseRoot original root found by a previous non-bracketing solver
93       * @param min minimal bound of the search interval
94       * @param max maximal bound of the search interval
95       * @param allowedSolution the kind of solutions that the root-finding algorithm may
96       * accept as solutions.
97       * @return a root approximation, on the specified side of the exact root
98       * @throws NoBracketingException if the function has the same sign at the
99       * endpoints.
100      */
101     public static double forceSide(final int maxEval, final UnivariateFunction f,
102                                    final BracketedUnivariateSolver<UnivariateFunction> bracketing,
103                                    final double baseRoot, final double min, final double max,
104                                    final AllowedSolution allowedSolution)
105         throws NoBracketingException {
106 
107         if (allowedSolution == AllowedSolution.ANY_SIDE) {
108             // no further bracketing required
109             return baseRoot;
110         }
111 
112         // find a very small interval bracketing the root
113         final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
114                                          FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
115         double xLo        = FastMath.max(min, baseRoot - step);
116         double fLo        = f.value(xLo);
117         double xHi        = FastMath.min(max, baseRoot + step);
118         double fHi        = f.value(xHi);
119         int remainingEval = maxEval - 2;
120         while (remainingEval > 0) {
121 
122             if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
123                 // compute the root on the selected side
124                 return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
125             }
126 
127             // try increasing the interval
128             boolean changeLo = false;
129             boolean changeHi = false;
130             if (fLo < fHi) {
131                 // increasing function
132                 if (fLo >= 0) {
133                     changeLo = true;
134                 } else {
135                     changeHi = true;
136                 }
137             } else if (fLo > fHi) {
138                 // decreasing function
139                 if (fLo <= 0) {
140                     changeLo = true;
141                 } else {
142                     changeHi = true;
143                 }
144             } else {
145                 // unknown variation
146                 changeLo = true;
147                 changeHi = true;
148             }
149 
150             // update the lower bound
151             if (changeLo) {
152                 xLo = FastMath.max(min, xLo - step);
153                 fLo  = f.value(xLo);
154                 remainingEval--;
155             }
156 
157             // update the higher bound
158             if (changeHi) {
159                 xHi = FastMath.min(max, xHi + step);
160                 fHi  = f.value(xHi);
161                 remainingEval--;
162             }
163 
164         }
165 
166         throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
167                                         xLo, xHi, fLo, fHi,
168                                         maxEval - remainingEval, maxEval, baseRoot,
169                                         min, max);
170 
171     }
172 
173     /**
174      * This method attempts to find two values a and b satisfying <ul>
175      * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
176      * <li> <code> f(a) * f(b) < 0 </code></li>
177      * </ul>
178      * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
179      * and <code>b</code> bracket a root of f.
180      * <p>
181      * The algorithm starts by setting
182      * <code>a := initial -1; b := initial +1,</code> examines the value of the
183      * function at <code>a</code> and <code>b</code> and keeps moving
184      * the endpoints out by one unit each time through a loop that terminates
185      * when one of the following happens: <ul>
186      * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
187      * <li> <code> a = lower </code> and <code> b = upper</code>
188      * -- NoBracketingException </li>
189      * <li> <code> Integer.MAX_VALUE</code> iterations elapse
190      * -- NoBracketingException </li>
191      * </ul></p>
192      * <p>
193      * <strong>Note: </strong> this method can take
194      * <code>Integer.MAX_VALUE</code> iterations to throw a
195      * <code>ConvergenceException.</code>  Unless you are confident that there
196      * is a root between <code>lowerBound</code> and <code>upperBound</code>
197      * near <code>initial,</code> it is better to use
198      * {@link #bracket(UnivariateFunction, double, double, double, int)},
199      * explicitly specifying the maximum number of iterations.</p>
200      *
201      * @param function Function.
202      * @param initial Initial midpoint of interval being expanded to
203      * bracket a root.
204      * @param lowerBound Lower bound (a is never lower than this value)
205      * @param upperBound Upper bound (b never is greater than this
206      * value).
207      * @return a two-element array holding a and b.
208      * @throws NoBracketingException if a root cannot be bracketted.
209      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
210      * @throws NullArgumentException if {@code function} is {@code null}.
211      */
212     public static double[] bracket(UnivariateFunction function,
213                                    double initial,
214                                    double lowerBound, double upperBound)
215         throws NullArgumentException,
216                NotStrictlyPositiveException,
217                NoBracketingException {
218         return bracket(function, initial, lowerBound, upperBound, Integer.MAX_VALUE);
219     }
220 
221      /**
222      * This method attempts to find two values a and b satisfying <ul>
223      * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
224      * <li> <code> f(a) * f(b) <= 0 </code> </li>
225      * </ul>
226      * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
227      * and <code>b</code> bracket a root of f.
228      * <p>
229      * The algorithm starts by setting
230      * <code>a := initial -1; b := initial +1,</code> examines the value of the
231      * function at <code>a</code> and <code>b</code> and keeps moving
232      * the endpoints out by one unit each time through a loop that terminates
233      * when one of the following happens: <ul>
234      * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
235      * <li> <code> a = lower </code> and <code> b = upper</code>
236      * -- NoBracketingException </li>
237      * <li> <code> maximumIterations</code> iterations elapse
238      * -- NoBracketingException </li></ul></p>
239      *
240      * @param function Function.
241      * @param initial Initial midpoint of interval being expanded to
242      * bracket a root.
243      * @param lowerBound Lower bound (a is never lower than this value).
244      * @param upperBound Upper bound (b never is greater than this
245      * value).
246      * @param maximumIterations Maximum number of iterations to perform
247      * @return a two element array holding a and b.
248      * @throws NoBracketingException if the algorithm fails to find a and b
249      * satisfying the desired conditions.
250      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
251      * @throws NullArgumentException if {@code function} is {@code null}.
252      */
253     public static double[] bracket(UnivariateFunction function,
254                                    double initial,
255                                    double lowerBound, double upperBound,
256                                    int maximumIterations)
257         throws NullArgumentException,
258                NotStrictlyPositiveException,
259                NoBracketingException {
260         if (function == null) {
261             throw new NullArgumentException(LocalizedFormats.FUNCTION);
262         }
263         if (maximumIterations <= 0)  {
264             throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
265         }
266         verifySequence(lowerBound, initial, upperBound);
267 
268         double a = initial;
269         double b = initial;
270         double fa;
271         double fb;
272         int numIterations = 0;
273 
274         do {
275             a = FastMath.max(a - 1.0, lowerBound);
276             b = FastMath.min(b + 1.0, upperBound);
277             fa = function.value(a);
278 
279             fb = function.value(b);
280             ++numIterations;
281         } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
282                 ((a > lowerBound) || (b < upperBound)));
283 
284         if (fa * fb > 0.0) {
285             throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
286                                             a, b, fa, fb,
287                                             numIterations, maximumIterations, initial,
288                                             lowerBound, upperBound);
289         }
290 
291         return new double[] {a, b};
292     }
293 
294     /**
295      * Compute the midpoint of two values.
296      *
297      * @param a first value.
298      * @param b second value.
299      * @return the midpoint.
300      */
301     public static double midpoint(double a, double b) {
302         return (a + b) * 0.5;
303     }
304 
305     /**
306      * Check whether the interval bounds bracket a root. That is, if the
307      * values at the endpoints are not equal to zero, then the function takes
308      * opposite signs at the endpoints.
309      *
310      * @param function Function.
311      * @param lower Lower endpoint.
312      * @param upper Upper endpoint.
313      * @return {@code true} if the function values have opposite signs at the
314      * given points.
315      * @throws NullArgumentException if {@code function} is {@code null}.
316      */
317     public static boolean isBracketing(UnivariateFunction function,
318                                        final double lower,
319                                        final double upper)
320         throws NullArgumentException {
321         if (function == null) {
322             throw new NullArgumentException(LocalizedFormats.FUNCTION);
323         }
324         final double fLo = function.value(lower);
325         final double fHi = function.value(upper);
326         return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
327     }
328 
329     /**
330      * Check whether the arguments form a (strictly) increasing sequence.
331      *
332      * @param start First number.
333      * @param mid Second number.
334      * @param end Third number.
335      * @return {@code true} if the arguments form an increasing sequence.
336      */
337     public static boolean isSequence(final double start,
338                                      final double mid,
339                                      final double end) {
340         return (start < mid) && (mid < end);
341     }
342 
343     /**
344      * Check that the endpoints specify an interval.
345      *
346      * @param lower Lower endpoint.
347      * @param upper Upper endpoint.
348      * @throws NumberIsTooLargeException if {@code lower >= upper}.
349      */
350     public static void verifyInterval(final double lower,
351                                       final double upper)
352         throws NumberIsTooLargeException {
353         if (lower >= upper) {
354             throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
355                                                 lower, upper, false);
356         }
357     }
358 
359     /**
360      * Check that {@code lower < initial < upper}.
361      *
362      * @param lower Lower endpoint.
363      * @param initial Initial value.
364      * @param upper Upper endpoint.
365      * @throws NumberIsTooLargeException if {@code lower >= initial} or
366      * {@code initial >= upper}.
367      */
368     public static void verifySequence(final double lower,
369                                       final double initial,
370                                       final double upper)
371         throws NumberIsTooLargeException {
372         verifyInterval(lower, initial);
373         verifyInterval(initial, upper);
374     }
375 
376     /**
377      * Check that the endpoints specify an interval and the end points
378      * bracket a root.
379      *
380      * @param function Function.
381      * @param lower Lower endpoint.
382      * @param upper Upper endpoint.
383      * @throws NoBracketingException if the function has the same sign at the
384      * endpoints.
385      * @throws NullArgumentException if {@code function} is {@code null}.
386      */
387     public static void verifyBracketing(UnivariateFunction function,
388                                         final double lower,
389                                         final double upper)
390         throws NullArgumentException,
391                NoBracketingException {
392         if (function == null) {
393             throw new NullArgumentException(LocalizedFormats.FUNCTION);
394         }
395         verifyInterval(lower, upper);
396         if (!isBracketing(function, lower, upper)) {
397             throw new NoBracketingException(lower, upper,
398                                             function.value(lower),
399                                             function.value(upper));
400         }
401     }
402 }