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 }