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 1579346 2014-03-19 18:43:39Z 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 simply calls {@link #bracket(UnivariateFunction, double, double, double,
175      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
176      * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
177      * <strong>Note: </strong> this method can take
178      * <code>Integer.MAX_VALUE</code> iterations to throw a
179      * <code>ConvergenceException.</code>  Unless you are confident that there
180      * is a root between <code>lowerBound</code> and <code>upperBound</code>
181      * near <code>initial,</code> it is better to use
182      * {@link #bracket(UnivariateFunction, double, double, double, double,
183      * double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
184      * explicitly specifying the maximum number of iterations.</p>
185      *
186      * @param function Function.
187      * @param initial Initial midpoint of interval being expanded to
188      * bracket a root.
189      * @param lowerBound Lower bound (a is never lower than this value)
190      * @param upperBound Upper bound (b never is greater than this
191      * value).
192      * @return a two-element array holding a and b.
193      * @throws NoBracketingException if a root cannot be bracketted.
194      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
195      * @throws NullArgumentException if {@code function} is {@code null}.
196      */
197     public static double[] bracket(UnivariateFunction function,
198                                    double initial,
199                                    double lowerBound, double upperBound)
200         throws NullArgumentException,
201                NotStrictlyPositiveException,
202                NoBracketingException {
203         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, Integer.MAX_VALUE);
204     }
205 
206      /**
207      * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
208      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
209      * with {@code q} and {@code r} set to 1.0.
210      * @param function Function.
211      * @param initial Initial midpoint of interval being expanded to
212      * bracket a root.
213      * @param lowerBound Lower bound (a is never lower than this value).
214      * @param upperBound Upper bound (b never is greater than this
215      * value).
216      * @param maximumIterations Maximum number of iterations to perform
217      * @return a two element array holding a and b.
218      * @throws NoBracketingException if the algorithm fails to find a and b
219      * satisfying the desired conditions.
220      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
221      * @throws NullArgumentException if {@code function} is {@code null}.
222      */
223     public static double[] bracket(UnivariateFunction function,
224                                    double initial,
225                                    double lowerBound, double upperBound,
226                                    int maximumIterations)
227         throws NullArgumentException,
228                NotStrictlyPositiveException,
229                NoBracketingException {
230         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations);
231     }
232 
233     /**
234      * This method attempts to find two values a and b satisfying <ul>
235      * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
236      * <li> {@code f(a) * f(b) <= 0} </li>
237      * </ul>
238      * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
239      * and {@code b} bracket a root of {@code f}.
240      * <p>
241      * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
242      * values of k, where \( l_k = max(lower, initial - \delta_k) \),
243      * \( u_k = min(upper, initial + \delta_k) \), using recurrence
244      * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
245      * The algorithm stops when one of the following happens: <ul>
246      * <li> at least one positive and one negative value have been found --  success!</li>
247      * <li> both endpoints have reached their respective limites -- NoBracketingException </li>
248      * <li> {@code maximumIterations} iterations elapse -- NoBracketingException </li></ul></p>
249      * <p>
250      * If different signs are found at first iteration ({@code k=1}), then the returned
251      * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
252      * iteration ({code k>1}, then the returned interval will be either
253      * \( [a, b] = [l_{k+1}, l_{k}] \) or ( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
254      * with these parameters will therefore start with the smallest bracketing interval known
255      * at this step.
256      * </p>
257      * <p>
258      * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
259      * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
260      * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
261      * is larger than 1, the sequence has an asymtotically exponential rate. Note than the
262      * additive parameter {@code q} should never be set to zero, otherwise the interval would
263      * degenerate to the single initial point for all values of {@code k}.
264      * </p>
265      * <p>
266      * As a rule of thumb, when the location of the root is expected to be approximately known
267      * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
268      * order of magnitude of the error margin. When the location of the root is really a wild guess,
269      * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
270      * length at each iteration) and {@code q} should be set according to half the initial
271      * search interval length.
272      * </p>
273      * <p>
274      * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
275      * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
276      * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
277      * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
278      * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
279      * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
280      * bracketing interval.
281      * </p>
282      * @param function function to check
283      * @param initial Initial midpoint of interval being expanded to
284      * bracket a root.
285      * @param lowerBound Lower bound (a is never lower than this value).
286      * @param upperBound Upper bound (b never is greater than this
287      * value).
288      * @param q additive offset used to compute bounds sequence (must be strictly positive)
289      * @param r multiplicative factor used to compute bounds sequence
290      * @param maximumIterations Maximum number of iterations to perform
291      * @return a two element array holding the bracketing values.
292      * @exception NoBracketingException if function cannot be bracketed in the search interval
293      */
294     public static double[] bracket(final UnivariateFunction function, final double initial,
295                                    final double lowerBound, final double upperBound,
296                                    final double q, final double r, final int maximumIterations)
297         throws NoBracketingException {
298 
299         if (function == null) {
300             throw new NullArgumentException(LocalizedFormats.FUNCTION);
301         }
302         if (q <= 0)  {
303             throw new NotStrictlyPositiveException(q);
304         }
305         if (maximumIterations <= 0)  {
306             throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
307         }
308         verifySequence(lowerBound, initial, upperBound);
309 
310         // initialize the recurrence
311         double a     = initial;
312         double b     = initial;
313         double fa    = Double.NaN;
314         double fb    = Double.NaN;
315         double delta = 0;
316 
317         for (int numIterations = 0;
318              (numIterations < maximumIterations) && (a > lowerBound || b > upperBound);
319              ++numIterations) {
320 
321             final double previousA  = a;
322             final double previousFa = fa;
323             final double previousB  = b;
324             final double previousFb = fb;
325 
326             delta = r * delta + q;
327             a     = FastMath.max(initial - delta, lowerBound);
328             b     = FastMath.min(initial + delta, upperBound);
329             fa    = function.value(a);
330             fb    = function.value(b);
331 
332             if (numIterations == 0) {
333                 // at first iteration, we don't have a previous interval
334                 // we simply compare both sides of the initial interval
335                 if (fa * fb <= 0) {
336                     // the first interval already brackets a root
337                     return new double[] { a, b };
338                 }
339             } else {
340                 // we have a previous interval with constant sign and expand it,
341                 // we expect sign changes to occur at boundaries
342                 if (fa * previousFa <= 0) {
343                     // sign change detected at near lower bound
344                     return new double[] { a, previousA };
345                 } else if (fb * previousFb <= 0) {
346                     // sign change detected at near upper bound
347                     return new double[] { previousB, b };
348                 }
349             }
350 
351         }
352 
353         // no bracketing found
354         throw new NoBracketingException(a, b, fa, fb);
355 
356     }
357 
358     /**
359      * Compute the midpoint of two values.
360      *
361      * @param a first value.
362      * @param b second value.
363      * @return the midpoint.
364      */
365     public static double midpoint(double a, double b) {
366         return (a + b) * 0.5;
367     }
368 
369     /**
370      * Check whether the interval bounds bracket a root. That is, if the
371      * values at the endpoints are not equal to zero, then the function takes
372      * opposite signs at the endpoints.
373      *
374      * @param function Function.
375      * @param lower Lower endpoint.
376      * @param upper Upper endpoint.
377      * @return {@code true} if the function values have opposite signs at the
378      * given points.
379      * @throws NullArgumentException if {@code function} is {@code null}.
380      */
381     public static boolean isBracketing(UnivariateFunction function,
382                                        final double lower,
383                                        final double upper)
384         throws NullArgumentException {
385         if (function == null) {
386             throw new NullArgumentException(LocalizedFormats.FUNCTION);
387         }
388         final double fLo = function.value(lower);
389         final double fHi = function.value(upper);
390         return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
391     }
392 
393     /**
394      * Check whether the arguments form a (strictly) increasing sequence.
395      *
396      * @param start First number.
397      * @param mid Second number.
398      * @param end Third number.
399      * @return {@code true} if the arguments form an increasing sequence.
400      */
401     public static boolean isSequence(final double start,
402                                      final double mid,
403                                      final double end) {
404         return (start < mid) && (mid < end);
405     }
406 
407     /**
408      * Check that the endpoints specify an interval.
409      *
410      * @param lower Lower endpoint.
411      * @param upper Upper endpoint.
412      * @throws NumberIsTooLargeException if {@code lower >= upper}.
413      */
414     public static void verifyInterval(final double lower,
415                                       final double upper)
416         throws NumberIsTooLargeException {
417         if (lower >= upper) {
418             throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
419                                                 lower, upper, false);
420         }
421     }
422 
423     /**
424      * Check that {@code lower < initial < upper}.
425      *
426      * @param lower Lower endpoint.
427      * @param initial Initial value.
428      * @param upper Upper endpoint.
429      * @throws NumberIsTooLargeException if {@code lower >= initial} or
430      * {@code initial >= upper}.
431      */
432     public static void verifySequence(final double lower,
433                                       final double initial,
434                                       final double upper)
435         throws NumberIsTooLargeException {
436         verifyInterval(lower, initial);
437         verifyInterval(initial, upper);
438     }
439 
440     /**
441      * Check that the endpoints specify an interval and the end points
442      * bracket a root.
443      *
444      * @param function Function.
445      * @param lower Lower endpoint.
446      * @param upper Upper endpoint.
447      * @throws NoBracketingException if the function has the same sign at the
448      * endpoints.
449      * @throws NullArgumentException if {@code function} is {@code null}.
450      */
451     public static void verifyBracketing(UnivariateFunction function,
452                                         final double lower,
453                                         final double upper)
454         throws NullArgumentException,
455                NoBracketingException {
456         if (function == null) {
457             throw new NullArgumentException(LocalizedFormats.FUNCTION);
458         }
459         verifyInterval(lower, upper);
460         if (!isBracketing(function, lower, upper)) {
461             throw new NoBracketingException(lower, upper,
462                                             function.value(lower),
463                                             function.value(upper));
464         }
465     }
466 }