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.math3.analysis.solvers;
018    
019    import org.apache.commons.math3.analysis.UnivariateFunction;
020    import org.apache.commons.math3.exception.NoBracketingException;
021    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
022    import org.apache.commons.math3.exception.NullArgumentException;
023    import org.apache.commons.math3.exception.NumberIsTooLargeException;
024    import org.apache.commons.math3.exception.util.LocalizedFormats;
025    import org.apache.commons.math3.util.FastMath;
026    
027    /**
028     * Utility routines for {@link UnivariateSolver} objects.
029     *
030     * @version $Id: UnivariateSolverUtils.java 1400850 2012-10-22 11:57:17Z erans $
031     */
032    public class UnivariateSolverUtils {
033        /**
034         * Class contains only static methods.
035         */
036        private UnivariateSolverUtils() {}
037    
038        /**
039         * Convenience method to find a zero of a univariate real function.  A default
040         * solver is used.
041         *
042         * @param function Function.
043         * @param x0 Lower bound for the interval.
044         * @param x1 Upper bound for the interval.
045         * @return a value where the function is zero.
046         * @throws NoBracketingException if the function has the same sign at the
047         * endpoints.
048         * @throws NullArgumentException if {@code function} is {@code null}.
049         */
050        public static double solve(UnivariateFunction function, double x0, double x1)
051            throws NullArgumentException,
052                   NoBracketingException {
053            if (function == null) {
054                throw new NullArgumentException(LocalizedFormats.FUNCTION);
055            }
056            final UnivariateSolver solver = new BrentSolver();
057            return solver.solve(Integer.MAX_VALUE, function, x0, x1);
058        }
059    
060        /**
061         * Convenience method to find a zero of a univariate real function.  A default
062         * solver is used.
063         *
064         * @param function Function.
065         * @param x0 Lower bound for the interval.
066         * @param x1 Upper bound for the interval.
067         * @param absoluteAccuracy Accuracy to be used by the solver.
068         * @return a value where the function is zero.
069         * @throws NoBracketingException if the function has the same sign at the
070         * endpoints.
071         * @throws NullArgumentException if {@code function} is {@code null}.
072         */
073        public static double solve(UnivariateFunction function,
074                                   double x0, double x1,
075                                   double absoluteAccuracy)
076            throws NullArgumentException,
077                   NoBracketingException {
078            if (function == null) {
079                throw new NullArgumentException(LocalizedFormats.FUNCTION);
080            }
081            final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
082            return solver.solve(Integer.MAX_VALUE, function, x0, x1);
083        }
084    
085        /** Force a root found by a non-bracketing solver to lie on a specified side,
086         * as if the solver was a bracketing one.
087         * @param maxEval maximal number of new evaluations of the function
088         * (evaluations already done for finding the root should have already been subtracted
089         * from this number)
090         * @param f function to solve
091         * @param bracketing bracketing solver to use for shifting the root
092         * @param baseRoot original root found by a previous non-bracketing solver
093         * @param min minimal bound of the search interval
094         * @param max maximal bound of the search interval
095         * @param allowedSolution the kind of solutions that the root-finding algorithm may
096         * accept as solutions.
097         * @return a root approximation, on the specified side of the exact root
098         * @throws NoBracketingException if the function has the same sign at the
099         * 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    }