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   */
31  public class UnivariateSolverUtils {
32      /**
33       * Class contains only static methods.
34       */
35      private UnivariateSolverUtils() {}
36  
37      /**
38       * Convenience method to find a zero of a univariate real function.  A default
39       * solver is used.
40       *
41       * @param function Function.
42       * @param x0 Lower bound for the interval.
43       * @param x1 Upper bound for the interval.
44       * @return a value where the function is zero.
45       * @throws NoBracketingException if the function has the same sign at the
46       * endpoints.
47       * @throws NullArgumentException if {@code function} is {@code null}.
48       */
49      public static double solve(UnivariateFunction function, double x0, double x1)
50          throws NullArgumentException,
51                 NoBracketingException {
52          if (function == null) {
53              throw new NullArgumentException(LocalizedFormats.FUNCTION);
54          }
55          final UnivariateSolver solver = new BrentSolver();
56          return solver.solve(Integer.MAX_VALUE, function, x0, x1);
57      }
58  
59      /**
60       * Convenience method to find a zero of a univariate real function.  A default
61       * solver is used.
62       *
63       * @param function Function.
64       * @param x0 Lower bound for the interval.
65       * @param x1 Upper bound for the interval.
66       * @param absoluteAccuracy Accuracy to be used by the solver.
67       * @return a value where the function is zero.
68       * @throws NoBracketingException if the function has the same sign at the
69       * endpoints.
70       * @throws NullArgumentException if {@code function} is {@code null}.
71       */
72      public static double solve(UnivariateFunction function,
73                                 double x0, double x1,
74                                 double absoluteAccuracy)
75          throws NullArgumentException,
76                 NoBracketingException {
77          if (function == null) {
78              throw new NullArgumentException(LocalizedFormats.FUNCTION);
79          }
80          final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
81          return solver.solve(Integer.MAX_VALUE, function, x0, x1);
82      }
83  
84      /** Force a root found by a non-bracketing solver to lie on a specified side,
85       * as if the solver was a bracketing one.
86       * @param maxEval maximal number of new evaluations of the function
87       * (evaluations already done for finding the root should have already been subtracted
88       * from this number)
89       * @param f function to solve
90       * @param bracketing bracketing solver to use for shifting the root
91       * @param baseRoot original root found by a previous non-bracketing solver
92       * @param min minimal bound of the search interval
93       * @param max maximal bound of the search interval
94       * @param allowedSolution the kind of solutions that the root-finding algorithm may
95       * accept as solutions.
96       * @return a root approximation, on the specified side of the exact root
97       * @throws NoBracketingException if the function has the same sign at the
98       * endpoints.
99       */
100     public static double forceSide(final int maxEval, final UnivariateFunction f,
101                                    final BracketedUnivariateSolver<UnivariateFunction> bracketing,
102                                    final double baseRoot, final double min, final double max,
103                                    final AllowedSolution allowedSolution)
104         throws NoBracketingException {
105 
106         if (allowedSolution == AllowedSolution.ANY_SIDE) {
107             // no further bracketing required
108             return baseRoot;
109         }
110 
111         // find a very small interval bracketing the root
112         final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
113                                          FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
114         double xLo        = FastMath.max(min, baseRoot - step);
115         double fLo        = f.value(xLo);
116         double xHi        = FastMath.min(max, baseRoot + step);
117         double fHi        = f.value(xHi);
118         int remainingEval = maxEval - 2;
119         while (remainingEval > 0) {
120 
121             if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
122                 // compute the root on the selected side
123                 return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
124             }
125 
126             // try increasing the interval
127             boolean changeLo = false;
128             boolean changeHi = false;
129             if (fLo < fHi) {
130                 // increasing function
131                 if (fLo >= 0) {
132                     changeLo = true;
133                 } else {
134                     changeHi = true;
135                 }
136             } else if (fLo > fHi) {
137                 // decreasing function
138                 if (fLo <= 0) {
139                     changeLo = true;
140                 } else {
141                     changeHi = true;
142                 }
143             } else {
144                 // unknown variation
145                 changeLo = true;
146                 changeHi = true;
147             }
148 
149             // update the lower bound
150             if (changeLo) {
151                 xLo = FastMath.max(min, xLo - step);
152                 fLo  = f.value(xLo);
153                 remainingEval--;
154             }
155 
156             // update the higher bound
157             if (changeHi) {
158                 xHi = FastMath.min(max, xHi + step);
159                 fHi  = f.value(xHi);
160                 remainingEval--;
161             }
162 
163         }
164 
165         throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
166                                         xLo, xHi, fLo, fHi,
167                                         maxEval - remainingEval, maxEval, baseRoot,
168                                         min, max);
169 
170     }
171 
172     /**
173      * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
174      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
175      * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
176      * <strong>Note: </strong> this method can take
177      * <code>Integer.MAX_VALUE</code> iterations to throw a
178      * <code>ConvergenceException.</code>  Unless you are confident that there
179      * is a root between <code>lowerBound</code> and <code>upperBound</code>
180      * near <code>initial,</code> it is better to use
181      * {@link #bracket(UnivariateFunction, double, double, double, double,
182      * double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
183      * explicitly specifying the maximum number of iterations.</p>
184      *
185      * @param function Function.
186      * @param initial Initial midpoint of interval being expanded to
187      * bracket a root.
188      * @param lowerBound Lower bound (a is never lower than this value)
189      * @param upperBound Upper bound (b never is greater than this
190      * value).
191      * @return a two-element array holding a and b.
192      * @throws NoBracketingException if a root cannot be bracketted.
193      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
194      * @throws NullArgumentException if {@code function} is {@code null}.
195      */
196     public static double[] bracket(UnivariateFunction function,
197                                    double initial,
198                                    double lowerBound, double upperBound)
199         throws NullArgumentException,
200                NotStrictlyPositiveException,
201                NoBracketingException {
202         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, Integer.MAX_VALUE);
203     }
204 
205      /**
206      * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
207      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
208      * with {@code q} and {@code r} set to 1.0.
209      * @param function Function.
210      * @param initial Initial midpoint of interval being expanded to
211      * bracket a root.
212      * @param lowerBound Lower bound (a is never lower than this value).
213      * @param upperBound Upper bound (b never is greater than this
214      * value).
215      * @param maximumIterations Maximum number of iterations to perform
216      * @return a two element array holding a and b.
217      * @throws NoBracketingException if the algorithm fails to find a and b
218      * satisfying the desired conditions.
219      * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}.
220      * @throws NullArgumentException if {@code function} is {@code null}.
221      */
222     public static double[] bracket(UnivariateFunction function,
223                                    double initial,
224                                    double lowerBound, double upperBound,
225                                    int maximumIterations)
226         throws NullArgumentException,
227                NotStrictlyPositiveException,
228                NoBracketingException {
229         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations);
230     }
231 
232     /**
233      * This method attempts to find two values a and b satisfying <ul>
234      * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
235      * <li> {@code f(a) * f(b) <= 0} </li>
236      * </ul>
237      * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
238      * and {@code b} bracket a root of {@code f}.
239      * <p>
240      * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
241      * values of k, where \( l_k = max(lower, initial - \delta_k) \),
242      * \( u_k = min(upper, initial + \delta_k) \), using recurrence
243      * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
244      * The algorithm stops when one of the following happens: <ul>
245      * <li> at least one positive and one negative value have been found --  success!</li>
246      * <li> both endpoints have reached their respective limites -- NoBracketingException </li>
247      * <li> {@code maximumIterations} iterations elapse -- NoBracketingException </li></ul></p>
248      * <p>
249      * If different signs are found at first iteration ({@code k=1}), then the returned
250      * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
251      * iteration ({code k>1}, then the returned interval will be either
252      * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
253      * with these parameters will therefore start with the smallest bracketing interval known
254      * at this step.
255      * </p>
256      * <p>
257      * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
258      * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
259      * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
260      * is larger than 1, the sequence has an asymtotically exponential rate. Note than the
261      * additive parameter {@code q} should never be set to zero, otherwise the interval would
262      * degenerate to the single initial point for all values of {@code k}.
263      * </p>
264      * <p>
265      * As a rule of thumb, when the location of the root is expected to be approximately known
266      * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
267      * order of magnitude of the error margin. When the location of the root is really a wild guess,
268      * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
269      * length at each iteration) and {@code q} should be set according to half the initial
270      * search interval length.
271      * </p>
272      * <p>
273      * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
274      * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
275      * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
276      * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
277      * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
278      * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
279      * bracketing interval.
280      * </p>
281      * @param function function to check
282      * @param initial Initial midpoint of interval being expanded to
283      * bracket a root.
284      * @param lowerBound Lower bound (a is never lower than this value).
285      * @param upperBound Upper bound (b never is greater than this
286      * value).
287      * @param q additive offset used to compute bounds sequence (must be strictly positive)
288      * @param r multiplicative factor used to compute bounds sequence
289      * @param maximumIterations Maximum number of iterations to perform
290      * @return a two element array holding the bracketing values.
291      * @exception NoBracketingException if function cannot be bracketed in the search interval
292      */
293     public static double[] bracket(final UnivariateFunction function, final double initial,
294                                    final double lowerBound, final double upperBound,
295                                    final double q, final double r, final int maximumIterations)
296         throws NoBracketingException {
297 
298         if (function == null) {
299             throw new NullArgumentException(LocalizedFormats.FUNCTION);
300         }
301         if (q <= 0)  {
302             throw new NotStrictlyPositiveException(q);
303         }
304         if (maximumIterations <= 0)  {
305             throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
306         }
307         verifySequence(lowerBound, initial, upperBound);
308 
309         // initialize the recurrence
310         double a     = initial;
311         double b     = initial;
312         double fa    = Double.NaN;
313         double fb    = Double.NaN;
314         double delta = 0;
315 
316         for (int numIterations = 0;
317              (numIterations < maximumIterations) && (a > lowerBound || b > upperBound);
318              ++numIterations) {
319 
320             final double previousA  = a;
321             final double previousFa = fa;
322             final double previousB  = b;
323             final double previousFb = fb;
324 
325             delta = r * delta + q;
326             a     = FastMath.max(initial - delta, lowerBound);
327             b     = FastMath.min(initial + delta, upperBound);
328             fa    = function.value(a);
329             fb    = function.value(b);
330 
331             if (numIterations == 0) {
332                 // at first iteration, we don't have a previous interval
333                 // we simply compare both sides of the initial interval
334                 if (fa * fb <= 0) {
335                     // the first interval already brackets a root
336                     return new double[] { a, b };
337                 }
338             } else {
339                 // we have a previous interval with constant sign and expand it,
340                 // we expect sign changes to occur at boundaries
341                 if (fa * previousFa <= 0) {
342                     // sign change detected at near lower bound
343                     return new double[] { a, previousA };
344                 } else if (fb * previousFb <= 0) {
345                     // sign change detected at near upper bound
346                     return new double[] { previousB, b };
347                 }
348             }
349 
350         }
351 
352         // no bracketing found
353         throw new NoBracketingException(a, b, fa, fb);
354 
355     }
356 
357     /**
358      * Compute the midpoint of two values.
359      *
360      * @param a first value.
361      * @param b second value.
362      * @return the midpoint.
363      */
364     public static double midpoint(double a, double b) {
365         return (a + b) * 0.5;
366     }
367 
368     /**
369      * Check whether the interval bounds bracket a root. That is, if the
370      * values at the endpoints are not equal to zero, then the function takes
371      * opposite signs at the endpoints.
372      *
373      * @param function Function.
374      * @param lower Lower endpoint.
375      * @param upper Upper endpoint.
376      * @return {@code true} if the function values have opposite signs at the
377      * given points.
378      * @throws NullArgumentException if {@code function} is {@code null}.
379      */
380     public static boolean isBracketing(UnivariateFunction function,
381                                        final double lower,
382                                        final double upper)
383         throws NullArgumentException {
384         if (function == null) {
385             throw new NullArgumentException(LocalizedFormats.FUNCTION);
386         }
387         final double fLo = function.value(lower);
388         final double fHi = function.value(upper);
389         return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
390     }
391 
392     /**
393      * Check whether the arguments form a (strictly) increasing sequence.
394      *
395      * @param start First number.
396      * @param mid Second number.
397      * @param end Third number.
398      * @return {@code true} if the arguments form an increasing sequence.
399      */
400     public static boolean isSequence(final double start,
401                                      final double mid,
402                                      final double end) {
403         return (start < mid) && (mid < end);
404     }
405 
406     /**
407      * Check that the endpoints specify an interval.
408      *
409      * @param lower Lower endpoint.
410      * @param upper Upper endpoint.
411      * @throws NumberIsTooLargeException if {@code lower >= upper}.
412      */
413     public static void verifyInterval(final double lower,
414                                       final double upper)
415         throws NumberIsTooLargeException {
416         if (lower >= upper) {
417             throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
418                                                 lower, upper, false);
419         }
420     }
421 
422     /**
423      * Check that {@code lower < initial < upper}.
424      *
425      * @param lower Lower endpoint.
426      * @param initial Initial value.
427      * @param upper Upper endpoint.
428      * @throws NumberIsTooLargeException if {@code lower >= initial} or
429      * {@code initial >= upper}.
430      */
431     public static void verifySequence(final double lower,
432                                       final double initial,
433                                       final double upper)
434         throws NumberIsTooLargeException {
435         verifyInterval(lower, initial);
436         verifyInterval(initial, upper);
437     }
438 
439     /**
440      * Check that the endpoints specify an interval and the end points
441      * bracket a root.
442      *
443      * @param function Function.
444      * @param lower Lower endpoint.
445      * @param upper Upper endpoint.
446      * @throws NoBracketingException if the function has the same sign at the
447      * endpoints.
448      * @throws NullArgumentException if {@code function} is {@code null}.
449      */
450     public static void verifyBracketing(UnivariateFunction function,
451                                         final double lower,
452                                         final double upper)
453         throws NullArgumentException,
454                NoBracketingException {
455         if (function == null) {
456             throw new NullArgumentException(LocalizedFormats.FUNCTION);
457         }
458         verifyInterval(lower, upper);
459         if (!isBracketing(function, lower, upper)) {
460             throw new NoBracketingException(lower, upper,
461                                             function.value(lower),
462                                             function.value(upper));
463         }
464     }
465 }