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 */
17package org.apache.commons.math.analysis.solvers;
18
19import org.apache.commons.math.ConvergenceException;
20import org.apache.commons.math.FunctionEvaluationException;
21import org.apache.commons.math.MathRuntimeException;
22import org.apache.commons.math.analysis.UnivariateRealFunction;
23import org.apache.commons.math.exception.util.LocalizedFormats;
24import org.apache.commons.math.exception.NullArgumentException;
25import org.apache.commons.math.util.FastMath;
26
27/**
28 * Utility routines for {@link UnivariateRealSolver} objects.
29 *
30 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
31 */
32public class UnivariateRealSolverUtils {
33
34    /**
35     * Default constructor.
36     */
37    private UnivariateRealSolverUtils() {
38        super();
39    }
40
41    /**
42     * Convenience method to find a zero of a univariate real function.  A default
43     * solver is used.
44     *
45     * @param f the function.
46     * @param x0 the lower bound for the interval.
47     * @param x1 the upper bound for the interval.
48     * @return a value where the function is zero.
49     * @throws ConvergenceException if the iteration count was exceeded
50     * @throws FunctionEvaluationException if an error occurs evaluating the function
51     * @throws IllegalArgumentException if f is null or the endpoints do not
52     * specify a valid interval
53     */
54    public static double solve(UnivariateRealFunction f, double x0, double x1)
55    throws ConvergenceException, FunctionEvaluationException {
56        setup(f);
57        return LazyHolder.FACTORY.newDefaultSolver().solve(f, 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 f the function
65     * @param x0 the lower bound for the interval
66     * @param x1 the upper bound for the interval
67     * @param absoluteAccuracy the accuracy to be used by the solver
68     * @return a value where the function is zero
69     * @throws ConvergenceException if the iteration count is exceeded
70     * @throws FunctionEvaluationException if an error occurs evaluating the function
71     * @throws IllegalArgumentException if f is null, the endpoints do not
72     * specify a valid interval, or the absoluteAccuracy is not valid for the
73     * default solver
74     */
75    public static double solve(UnivariateRealFunction f, double x0, double x1,
76            double absoluteAccuracy) throws ConvergenceException,
77            FunctionEvaluationException {
78
79        setup(f);
80        UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
81        solver.setAbsoluteAccuracy(absoluteAccuracy);
82        return solver.solve(f, x0, x1);
83    }
84
85    /**
86     * This method attempts to find two values a and b satisfying <ul>
87    * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
88     * <li> <code> f(a) * f(b) < 0 </code></li>
89     * </ul>
90     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
91     * and <code>b</code> bracket a root of f.
92     * <p>
93     * The algorithm starts by setting
94     * <code>a := initial -1; b := initial +1,</code> examines the value of the
95     * function at <code>a</code> and <code>b</code> and keeps moving
96     * the endpoints out by one unit each time through a loop that terminates
97     * when one of the following happens: <ul>
98     * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
99     * <li> <code> a = lower </code> and <code> b = upper</code>
100     * -- ConvergenceException </li>
101     * <li> <code> Integer.MAX_VALUE</code> iterations elapse
102     * -- ConvergenceException </li>
103     * </ul></p>
104     * <p>
105     * <strong>Note: </strong> this method can take
106     * <code>Integer.MAX_VALUE</code> iterations to throw a
107     * <code>ConvergenceException.</code>  Unless you are confident that there
108     * is a root between <code>lowerBound</code> and <code>upperBound</code>
109     * near <code>initial,</code> it is better to use
110     * {@link #bracket(UnivariateRealFunction, double, double, double, int)},
111     * explicitly specifying the maximum number of iterations.</p>
112     *
113     * @param function the function
114     * @param initial initial midpoint of interval being expanded to
115     * bracket a root
116     * @param lowerBound lower bound (a is never lower than this value)
117     * @param upperBound upper bound (b never is greater than this
118     * value)
119     * @return a two element array holding {a, b}
120     * @throws ConvergenceException if a root can not be bracketted
121     * @throws FunctionEvaluationException if an error occurs evaluating the function
122     * @throws IllegalArgumentException if function is null, maximumIterations
123     * is not positive, or initial is not between lowerBound and upperBound
124     */
125    public static double[] bracket(UnivariateRealFunction function,
126            double initial, double lowerBound, double upperBound)
127    throws ConvergenceException, FunctionEvaluationException {
128        return bracket( function, initial, lowerBound, upperBound,
129            Integer.MAX_VALUE ) ;
130    }
131
132     /**
133     * This method attempts to find two values a and b satisfying <ul>
134     * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
135     * <li> <code> f(a) * f(b) <= 0 </code> </li>
136     * </ul>
137     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
138     * and <code>b</code> bracket a root of f.
139     * <p>
140     * The algorithm starts by setting
141     * <code>a := initial -1; b := initial +1,</code> examines the value of the
142     * function at <code>a</code> and <code>b</code> and keeps moving
143     * the endpoints out by one unit each time through a loop that terminates
144     * when one of the following happens: <ul>
145     * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
146     * <li> <code> a = lower </code> and <code> b = upper</code>
147     * -- ConvergenceException </li>
148     * <li> <code> maximumIterations</code> iterations elapse
149     * -- ConvergenceException </li></ul></p>
150     *
151     * @param function the function
152     * @param initial initial midpoint of interval being expanded to
153     * bracket a root
154     * @param lowerBound lower bound (a is never lower than this value)
155     * @param upperBound upper bound (b never is greater than this
156     * value)
157     * @param maximumIterations maximum number of iterations to perform
158     * @return a two element array holding {a, b}.
159     * @throws ConvergenceException if the algorithm fails to find a and b
160     * satisfying the desired conditions
161     * @throws FunctionEvaluationException if an error occurs evaluating the 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 new NullArgumentException(LocalizedFormats.FUNCTION);
172        }
173        if (maximumIterations <= 0)  {
174            throw MathRuntimeException.createIllegalArgumentException(
175                  LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
176        }
177        if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
178            throw MathRuntimeException.createIllegalArgumentException(
179                  LocalizedFormats.INVALID_BRACKETING_PARAMETERS,
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 = FastMath.max(a - 1.0, lowerBound);
190            b = FastMath.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                      LocalizedFormats.FAILED_BRACKETING,
201                      numIterations, maximumIterations, initial,
202                      lowerBound, upperBound, a, b, fa, fb);
203        }
204
205        return new double[]{a, b};
206    }
207
208    /**
209     * Compute the midpoint of two values.
210     *
211     * @param a first value.
212     * @param b second value.
213     * @return the midpoint.
214     */
215    public static double midpoint(double a, double b) {
216        return (a + b) * .5;
217    }
218
219    /**
220     * Checks to see if f is null, throwing IllegalArgumentException if so.
221     * @param f  input function
222     * @throws IllegalArgumentException if f is null
223     */
224    private static void setup(UnivariateRealFunction f) {
225        if (f == null) {
226            throw new NullArgumentException(LocalizedFormats.FUNCTION);
227        }
228    }
229
230    // CHECKSTYLE: stop HideUtilityClassConstructor
231    /** Holder for the factory.
232     * <p>We use here the Initialization On Demand Holder Idiom.</p>
233     */
234    private static class LazyHolder {
235        /** Cached solver factory */
236        private static final UnivariateRealSolverFactory FACTORY = UnivariateRealSolverFactory.newInstance();
237    }
238    // CHECKSTYLE: resume HideUtilityClassConstructor
239
240}
241