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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.solvers;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.exception.LocalizedCoreFormats;
28 import org.hipparchus.exception.MathIllegalArgumentException;
29 import org.hipparchus.exception.NullArgumentException;
30 import org.hipparchus.util.FastMath;
31 import org.hipparchus.util.MathArrays;
32 import org.hipparchus.util.MathUtils;
33
34 /**
35 * Utility routines for {@link UnivariateSolver} objects.
36 *
37 */
38 public class UnivariateSolverUtils {
39 /**
40 * Class contains only static methods.
41 */
42 private UnivariateSolverUtils() {}
43
44 /**
45 * Convenience method to find a zero of a univariate real function. A default
46 * solver is used.
47 *
48 * @param function Function.
49 * @param x0 Lower bound for the interval.
50 * @param x1 Upper bound for the interval.
51 * @return a value where the function is zero.
52 * @throws MathIllegalArgumentException if the function has the same sign at the
53 * endpoints.
54 * @throws NullArgumentException if {@code function} is {@code null}.
55 */
56 public static double solve(UnivariateFunction function, double x0, double x1)
57 throws MathIllegalArgumentException, NullArgumentException {
58 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
59 final UnivariateSolver solver = new BrentSolver();
60 return solver.solve(Integer.MAX_VALUE, function, x0, x1);
61 }
62
63 /**
64 * Convenience method to find a zero of a univariate real function. A default
65 * solver is used.
66 *
67 * @param function Function.
68 * @param x0 Lower bound for the interval.
69 * @param x1 Upper bound for the interval.
70 * @param absoluteAccuracy Accuracy to be used by the solver.
71 * @return a value where the function is zero.
72 * @throws MathIllegalArgumentException if the function has the same sign at the
73 * endpoints.
74 * @throws NullArgumentException if {@code function} is {@code null}.
75 */
76 public static double solve(UnivariateFunction function,
77 double x0, double x1,
78 double absoluteAccuracy)
79 throws MathIllegalArgumentException, NullArgumentException {
80 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
81 final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
82 return solver.solve(Integer.MAX_VALUE, function, x0, x1);
83 }
84
85 /**
86 * Force a root found by a non-bracketing solver to lie on a specified side,
87 * as if the solver were a bracketing one.
88 *
89 * @param maxEval maximal number of new evaluations of the function
90 * (evaluations already done for finding the root should have already been subtracted
91 * from this number)
92 * @param f function to solve
93 * @param bracketing bracketing solver to use for shifting the root
94 * @param baseRoot original root found by a previous non-bracketing solver
95 * @param min minimal bound of the search interval
96 * @param max maximal bound of the search interval
97 * @param allowedSolution the kind of solutions that the root-finding algorithm may
98 * accept as solutions.
99 * @return a root approximation, on the specified side of the exact root
100 * @throws MathIllegalArgumentException if the function has the same sign at the
101 * endpoints.
102 */
103 public static double forceSide(final int maxEval, final UnivariateFunction f,
104 final BracketedUnivariateSolver<UnivariateFunction> bracketing,
105 final double baseRoot, final double min, final double max,
106 final AllowedSolution allowedSolution)
107 throws MathIllegalArgumentException {
108
109 if (allowedSolution == AllowedSolution.ANY_SIDE) {
110 // no further bracketing required
111 return baseRoot;
112 }
113
114 // find a very small interval bracketing the root
115 final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
116 FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
117 double xLo = FastMath.max(min, baseRoot - step);
118 double fLo = f.value(xLo);
119 double xHi = FastMath.min(max, baseRoot + step);
120 double fHi = f.value(xHi);
121 int remainingEval = maxEval - 2;
122 while (remainingEval > 0) {
123
124 if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
125 // compute the root on the selected side
126 return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
127 }
128
129 // try increasing the interval
130 boolean changeLo = false;
131 boolean changeHi = false;
132 if (fLo < fHi) {
133 // increasing function
134 if (fLo >= 0) {
135 changeLo = true;
136 } else {
137 changeHi = true;
138 }
139 } else if (fLo > fHi) {
140 // decreasing function
141 if (fLo <= 0) {
142 changeLo = true;
143 } else {
144 changeHi = true;
145 }
146 } else {
147 // unknown variation
148 changeLo = true;
149 changeHi = true;
150 }
151
152 // update the lower bound
153 if (changeLo) {
154 xLo = FastMath.max(min, xLo - step);
155 fLo = f.value(xLo);
156 remainingEval--;
157 }
158
159 // update the higher bound
160 if (changeHi) {
161 xHi = FastMath.min(max, xHi + step);
162 fHi = f.value(xHi);
163 remainingEval--;
164 }
165
166 }
167
168 throw new MathIllegalArgumentException(LocalizedCoreFormats.FAILED_BRACKETING,
169 xLo, xHi, fLo, fHi,
170 maxEval - remainingEval, maxEval, baseRoot,
171 min, max);
172
173 }
174
175 /**
176 * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
177 * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
178 * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
179 * <p>
180 * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE}
181 * iterations to throw a {@code MathIllegalStateException.} Unless you are
182 * confident that there is a root between {@code lowerBound} and
183 * {@code upperBound} near {@code initial}, it is better to use
184 * {@link #bracket(UnivariateFunction, double, double, double, double,double, int)
185 * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
186 * explicitly specifying the maximum number of iterations.</p>
187 *
188 * @param function Function.
189 * @param initial Initial midpoint of interval being expanded to
190 * bracket a root.
191 * @param lowerBound Lower bound (a is never lower than this value)
192 * @param upperBound Upper bound (b never is greater than this
193 * value).
194 * @return a two-element array holding a and b.
195 * @throws MathIllegalArgumentException if a root cannot be bracketed.
196 * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
197 * @throws NullArgumentException if {@code function} is {@code null}.
198 */
199 public static double[] bracket(UnivariateFunction function,
200 double initial,
201 double lowerBound, double upperBound)
202 throws MathIllegalArgumentException, NullArgumentException {
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 MathIllegalArgumentException if the algorithm fails to find a and b
219 * satisfying the desired conditions.
220 * @throws MathIllegalArgumentException 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 MathIllegalArgumentException, NullArgumentException {
228 return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations);
229 }
230
231 /**
232 * This method attempts to find two values a and b satisfying <ul>
233 * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
234 * <li> {@code f(a) * f(b) <= 0} </li>
235 * </ul>
236 * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
237 * and {@code b} bracket a root of {@code f}.
238 * <p>
239 * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
240 * values of k, where \( l_k = max(lower, initial - \delta_k) \),
241 * \( u_k = min(upper, initial + \delta_k) \), using recurrence
242 * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
243 * The algorithm stops when one of the following happens: <ul>
244 * <li> at least one positive and one negative value have been found -- success!</li>
245 * <li> both endpoints have reached their respective limits -- MathIllegalArgumentException </li>
246 * <li> {@code maximumIterations} iterations elapse -- MathIllegalArgumentException </li></ul>
247 * <p>
248 * If different signs are found at first iteration ({@code k=1}), then the returned
249 * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
250 * iteration {@code k>1}, then the returned interval will be either
251 * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
252 * with these parameters will therefore start with the smallest bracketing interval known
253 * at this step.
254 * </p>
255 * <p>
256 * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
257 * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
258 * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
259 * is larger than 1, the sequence has an asymptotically exponential rate. Note than the
260 * additive parameter {@code q} should never be set to zero, otherwise the interval would
261 * degenerate to the single initial point for all values of {@code k}.
262 * </p>
263 * <p>
264 * As a rule of thumb, when the location of the root is expected to be approximately known
265 * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
266 * order of magnitude of the error margin. When the location of the root is really a wild guess,
267 * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
268 * length at each iteration) and {@code q} should be set according to half the initial
269 * search interval length.
270 * </p>
271 * <p>
272 * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
273 * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
274 * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
275 * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
276 * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
277 * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
278 * bracketing interval.
279 * </p>
280 * @param function function to check
281 * @param initial Initial midpoint of interval being expanded to
282 * bracket a root.
283 * @param lowerBound Lower bound (a is never lower than this value).
284 * @param upperBound Upper bound (b never is greater than this
285 * value).
286 * @param q additive offset used to compute bounds sequence (must be strictly positive)
287 * @param r multiplicative factor used to compute bounds sequence
288 * @param maximumIterations Maximum number of iterations to perform
289 * @return a two element array holding the bracketing values.
290 * @exception MathIllegalArgumentException if function cannot be bracketed in the search interval
291 */
292 public static double[] bracket(final UnivariateFunction function, final double initial,
293 final double lowerBound, final double upperBound,
294 final double q, final double r, final int maximumIterations)
295 throws MathIllegalArgumentException {
296
297 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
298
299 if (q <= 0) {
300 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
301 q, 0);
302 }
303 if (maximumIterations <= 0) {
304 throw new MathIllegalArgumentException(LocalizedCoreFormats.INVALID_MAX_ITERATIONS, maximumIterations);
305 }
306 verifySequence(lowerBound, initial, upperBound);
307
308 // initialize the recurrence
309 double a = initial;
310 double b = initial;
311 double fa = Double.NaN;
312 double fb = Double.NaN;
313 double delta = 0;
314
315 for (int numIterations = 0;
316 (numIterations < maximumIterations) && (a > lowerBound || b < upperBound);
317 ++numIterations) {
318
319 final double previousA = a;
320 final double previousFa = fa;
321 final double previousB = b;
322 final double previousFb = fb;
323
324 delta = r * delta + q;
325 a = FastMath.max(initial - delta, lowerBound);
326 b = FastMath.min(initial + delta, upperBound);
327 fa = function.value(a);
328 fb = function.value(b);
329
330 if (numIterations == 0) {
331 // at first iteration, we don't have a previous interval
332 // we simply compare both sides of the initial interval
333 if (fa * fb <= 0) {
334 // the first interval already brackets a root
335 return new double[] { a, b };
336 }
337 } else {
338 // we have a previous interval with constant sign and expand it,
339 // we expect sign changes to occur at boundaries
340 if (fa * previousFa <= 0) {
341 // sign change detected at near lower bound
342 return new double[] { a, previousA };
343 } else if (fb * previousFb <= 0) {
344 // sign change detected at near upper bound
345 return new double[] { previousB, b };
346 }
347 }
348
349 }
350
351 // no bracketing found
352 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
353 a, b, fa, fb);
354
355 }
356
357 /**
358 * This method simply calls {@link #bracket(CalculusFieldUnivariateFunction,
359 * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement, CalculusFieldElement,
360 * CalculusFieldElement, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
361 * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
362 * <p>
363 * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE}
364 * iterations to throw a {@code MathIllegalStateException.} Unless you are
365 * confident that there is a root between {@code lowerBound} and
366 * {@code upperBound} near {@code initial}, it is better to use
367 * {@link #bracket(UnivariateFunction, double, double, double, double,double, int)
368 * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
369 * explicitly specifying the maximum number of iterations.</p>
370 *
371 * @param function Function.
372 * @param initial Initial midpoint of interval being expanded to
373 * bracket a root.
374 * @param lowerBound Lower bound (a is never lower than this value)
375 * @param upperBound Upper bound (b never is greater than this
376 * value).
377 * @param <T> type of the field elements
378 * @return a two-element array holding a and b.
379 * @throws MathIllegalArgumentException if a root cannot be bracketed.
380 * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
381 * @throws NullArgumentException if {@code function} is {@code null}.
382 * @since 1.2
383 */
384 public static <T extends CalculusFieldElement<T>> T[] bracket(CalculusFieldUnivariateFunction<T> function,
385 T initial,
386 T lowerBound, T upperBound)
387 throws MathIllegalArgumentException, NullArgumentException {
388 return bracket(function, initial, lowerBound, upperBound,
389 initial.getField().getOne(), initial.getField().getOne(),
390 Integer.MAX_VALUE);
391 }
392
393 /**
394 * This method simply calls {@link #bracket(CalculusFieldUnivariateFunction,
395 * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement, CalculusFieldElement,
396 * CalculusFieldElement, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
397 * with {@code q} and {@code r} set to 1.0.
398 * @param function Function.
399 * @param initial Initial midpoint of interval being expanded to
400 * bracket a root.
401 * @param lowerBound Lower bound (a is never lower than this value).
402 * @param upperBound Upper bound (b never is greater than this
403 * value).
404 * @param maximumIterations Maximum number of iterations to perform
405 * @param <T> type of the field elements
406 * @return a two element array holding a and b.
407 * @throws MathIllegalArgumentException if the algorithm fails to find a and b
408 * satisfying the desired conditions.
409 * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
410 * @throws NullArgumentException if {@code function} is {@code null}.
411 * @since 1.2
412 */
413 public static <T extends CalculusFieldElement<T>> T[] bracket(CalculusFieldUnivariateFunction<T> function,
414 T initial,
415 T lowerBound, T upperBound,
416 int maximumIterations)
417 throws MathIllegalArgumentException, NullArgumentException {
418 return bracket(function, initial, lowerBound, upperBound,
419 initial.getField().getOne(), initial.getField().getOne(),
420 maximumIterations);
421 }
422
423 /**
424 * This method attempts to find two values a and b satisfying <ul>
425 * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
426 * <li> {@code f(a) * f(b) <= 0} </li>
427 * </ul>
428 * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
429 * and {@code b} bracket a root of {@code f}.
430 * <p>
431 * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
432 * values of k, where \( l_k = max(lower, initial - \delta_k) \),
433 * \( u_k = min(upper, initial + \delta_k) \), using recurrence
434 * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
435 * The algorithm stops when one of the following happens: <ul>
436 * <li> at least one positive and one negative value have been found -- success!</li>
437 * <li> both endpoints have reached their respective limits -- MathIllegalArgumentException </li>
438 * <li> {@code maximumIterations} iterations elapse -- MathIllegalArgumentException </li></ul>
439 * <p>
440 * If different signs are found at first iteration ({@code k=1}), then the returned
441 * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
442 * iteration {@code k>1}, then the returned interval will be either
443 * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
444 * with these parameters will therefore start with the smallest bracketing interval known
445 * at this step.
446 * </p>
447 * <p>
448 * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
449 * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
450 * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
451 * is larger than 1, the sequence has an asymptotically exponential rate. Note than the
452 * additive parameter {@code q} should never be set to zero, otherwise the interval would
453 * degenerate to the single initial point for all values of {@code k}.
454 * </p>
455 * <p>
456 * As a rule of thumb, when the location of the root is expected to be approximately known
457 * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
458 * order of magnitude of the error margin. When the location of the root is really a wild guess,
459 * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
460 * length at each iteration) and {@code q} should be set according to half the initial
461 * search interval length.
462 * </p>
463 * <p>
464 * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
465 * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
466 * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
467 * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
468 * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
469 * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
470 * bracketing interval.
471 * </p>
472 * @param function function to check
473 * @param initial Initial midpoint of interval being expanded to
474 * bracket a root.
475 * @param lowerBound Lower bound (a is never lower than this value).
476 * @param upperBound Upper bound (b never is greater than this
477 * value).
478 * @param q additive offset used to compute bounds sequence (must be strictly positive)
479 * @param r multiplicative factor used to compute bounds sequence
480 * @param maximumIterations Maximum number of iterations to perform
481 * @param <T> type of the field elements
482 * @return a two element array holding the bracketing values.
483 * @exception MathIllegalArgumentException if function cannot be bracketed in the search interval
484 * @since 1.2
485 */
486 public static <T extends CalculusFieldElement<T>> T[] bracket(final CalculusFieldUnivariateFunction<T> function,
487 final T initial,
488 final T lowerBound, final T upperBound,
489 final T q, final T r,
490 final int maximumIterations)
491 throws MathIllegalArgumentException {
492
493 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
494
495 if (q.getReal() <= 0) {
496 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
497 q, 0);
498 }
499 if (maximumIterations <= 0) {
500 throw new MathIllegalArgumentException(LocalizedCoreFormats.INVALID_MAX_ITERATIONS, maximumIterations);
501 }
502 verifySequence(lowerBound.getReal(), initial.getReal(), upperBound.getReal());
503
504 // initialize the recurrence
505 T a = initial;
506 T b = initial;
507 T fa = null;
508 T fb = null;
509 T delta = initial.getField().getZero();
510
511 for (int numIterations = 0;
512 (numIterations < maximumIterations) &&
513 (a.getReal() > lowerBound.getReal() || b.getReal() < upperBound.getReal());
514 ++numIterations) {
515
516 final T previousA = a;
517 final T previousFa = fa;
518 final T previousB = b;
519 final T previousFb = fb;
520
521 delta = r.multiply(delta).add(q);
522 a = max(initial.subtract(delta), lowerBound);
523 b = min(initial.add(delta), upperBound);
524 fa = function.value(a);
525 fb = function.value(b);
526
527 if (numIterations == 0) {
528 // at first iteration, we don't have a previous interval
529 // we simply compare both sides of the initial interval
530 if (fa.multiply(fb).getReal() <= 0) {
531 // the first interval already brackets a root
532 final T[] interval = MathArrays.buildArray(initial.getField(), 2);
533 interval[0] = a;
534 interval[1] = b;
535 return interval;
536 }
537 } else {
538 // we have a previous interval with constant sign and expand it,
539 // we expect sign changes to occur at boundaries
540 if (fa.multiply(previousFa).getReal() <= 0) {
541 // sign change detected at near lower bound
542 final T[] interval = MathArrays.buildArray(initial.getField(), 2);
543 interval[0] = a;
544 interval[1] = previousA;
545 return interval;
546 } else if (fb.multiply(previousFb).getReal() <= 0) {
547 // sign change detected at near upper bound
548 final T[] interval = MathArrays.buildArray(initial.getField(), 2);
549 interval[0] = previousB;
550 interval[1] = b;
551 return interval;
552 }
553 }
554
555 }
556
557 // no bracketing found
558 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
559 a.getReal(), b.getReal(), fa.getReal(), fb.getReal());
560
561 }
562
563 /** Compute the maximum of two values
564 * @param a first value
565 * @param b second value
566 * @param <T> type of the field elements
567 * @return b if a is lesser or equal to b, a otherwise
568 * @since 1.2
569 */
570 private static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
571 return (a.subtract(b).getReal() <= 0) ? b : a;
572 }
573
574 /** Compute the minimum of two values
575 * @param a first value
576 * @param b second value
577 * @param <T> type of the field elements
578 * @return a if a is lesser or equal to b, b otherwise
579 * @since 1.2
580 */
581 private static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
582 return (a.subtract(b).getReal() <= 0) ? a : b;
583 }
584
585 /**
586 * Compute the midpoint of two values.
587 *
588 * @param a first value.
589 * @param b second value.
590 * @return the midpoint.
591 */
592 public static double midpoint(double a, double b) {
593 return (a + b) * 0.5;
594 }
595
596 /**
597 * Check whether the interval bounds bracket a root. That is, if the
598 * values at the endpoints are not equal to zero, then the function takes
599 * opposite signs at the endpoints.
600 *
601 * @param function Function.
602 * @param lower Lower endpoint.
603 * @param upper Upper endpoint.
604 * @return {@code true} if the function values have opposite signs at the
605 * given points.
606 * @throws NullArgumentException if {@code function} is {@code null}.
607 */
608 public static boolean isBracketing(UnivariateFunction function,
609 final double lower,
610 final double upper)
611 throws NullArgumentException {
612 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
613 final double fLo = function.value(lower);
614 final double fHi = function.value(upper);
615 return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
616 }
617
618 /**
619 * Check whether the arguments form a (strictly) increasing sequence.
620 *
621 * @param start First number.
622 * @param mid Second number.
623 * @param end Third number.
624 * @return {@code true} if the arguments form an increasing sequence.
625 */
626 public static boolean isSequence(final double start,
627 final double mid,
628 final double end) {
629 return (start < mid) && (mid < end);
630 }
631
632 /**
633 * Check that the endpoints specify an interval.
634 *
635 * @param lower Lower endpoint.
636 * @param upper Upper endpoint.
637 * @throws MathIllegalArgumentException if {@code lower >= upper}.
638 */
639 public static void verifyInterval(final double lower,
640 final double upper)
641 throws MathIllegalArgumentException {
642 if (lower >= upper) {
643 throw new MathIllegalArgumentException(LocalizedCoreFormats.ENDPOINTS_NOT_AN_INTERVAL,
644 lower, upper, false);
645 }
646 }
647
648 /**
649 * Check that {@code lower < initial < upper}.
650 *
651 * @param lower Lower endpoint.
652 * @param initial Initial value.
653 * @param upper Upper endpoint.
654 * @throws MathIllegalArgumentException if {@code lower >= initial} or
655 * {@code initial >= upper}.
656 */
657 public static void verifySequence(final double lower,
658 final double initial,
659 final double upper)
660 throws MathIllegalArgumentException {
661 verifyInterval(lower, initial);
662 verifyInterval(initial, upper);
663 }
664
665 /**
666 * Check that the endpoints specify an interval and the end points
667 * bracket a root.
668 *
669 * @param function Function.
670 * @param lower Lower endpoint.
671 * @param upper Upper endpoint.
672 * @throws MathIllegalArgumentException if the function has the same sign at the
673 * endpoints.
674 * @throws NullArgumentException if {@code function} is {@code null}.
675 */
676 public static void verifyBracketing(UnivariateFunction function,
677 final double lower,
678 final double upper)
679 throws MathIllegalArgumentException, NullArgumentException {
680 MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
681 verifyInterval(lower, upper);
682 if (!isBracketing(function, lower, upper)) {
683 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
684 lower, upper,
685 function.value(lower), function.value(upper));
686 }
687 }
688 }