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    *      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  
23  package org.hipparchus.analysis.solvers;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.exception.MathRuntimeException;
30  
31  /** Interface for {@link UnivariateSolver (univariate real) root-finding
32   * algorithms} that maintain a bracketed solution. There are several advantages
33   * to having such root-finding algorithms:
34   * <ul>
35   *  <li>The bracketed solution guarantees that the root is kept within the
36   *      interval. As such, these algorithms generally also guarantee
37   *      convergence.</li>
38   *  <li>The bracketed solution means that we have the opportunity to only
39   *      return roots that are greater than or equal to the actual root, or
40   *      are less than or equal to the actual root. That is, we can control
41   *      whether under-approximations and over-approximations are
42   *      {@link AllowedSolution allowed solutions}. Other root-finding
43   *      algorithms can usually only guarantee that the solution (the root that
44   *      was found) is around the actual root.</li>
45   * </ul>
46   *
47   * <p>For backwards compatibility, all root-finding algorithms must have
48   * {@link AllowedSolution#ANY_SIDE ANY_SIDE} as default for the allowed
49   * solutions.</p>
50   *
51   * @see AllowedSolution
52   * @param <T> the type of the field elements
53   */
54  public interface BracketedRealFieldUnivariateSolver<T extends CalculusFieldElement<T>> {
55  
56      /**
57       * Get the maximum number of function evaluations.
58       *
59       * @return the maximum number of function evaluations.
60       */
61      int getMaxEvaluations();
62  
63      /**
64       * Get the number of evaluations of the objective function.
65       * The number of evaluations corresponds to the last call to the
66       * {@code optimize} method. It is 0 if the method has not been
67       * called yet.
68       *
69       * @return the number of evaluations of the objective function.
70       */
71      int getEvaluations();
72  
73      /**
74       * Get the absolute accuracy of the solver.  Solutions returned by the
75       * solver should be accurate to this tolerance, i.e., if &epsilon; is the
76       * absolute accuracy of the solver and {@code v} is a value returned by
77       * one of the {@code solve} methods, then a root of the function should
78       * exist somewhere in the interval ({@code v} - &epsilon;, {@code v} + &epsilon;).
79       *
80       * @return the absolute accuracy.
81       */
82      T getAbsoluteAccuracy();
83  
84      /**
85       * Get the relative accuracy of the solver.  The contract for relative
86       * accuracy is the same as {@link #getAbsoluteAccuracy()}, but using
87       * relative, rather than absolute error.  If &rho; is the relative accuracy
88       * configured for a solver and {@code v} is a value returned, then a root
89       * of the function should exist somewhere in the interval
90       * ({@code v} - &rho; {@code v}, {@code v} + &rho; {@code v}).
91       *
92       * @return the relative accuracy.
93       */
94      T getRelativeAccuracy();
95  
96      /**
97       * Get the function value accuracy of the solver.  If {@code v} is
98       * a value returned by the solver for a function {@code f},
99       * then by contract, {@code |f(v)|} should be less than or equal to
100      * the function value accuracy configured for the solver.
101      *
102      * @return the function value accuracy.
103      */
104     T getFunctionValueAccuracy();
105 
106     /**
107      * Solve for a zero in the given interval.
108      * A solver may require that the interval brackets a single zero root.
109      * Solvers that do require bracketing should be able to handle the case
110      * where one of the endpoints is itself a root.
111      *
112      * @param maxEval Maximum number of evaluations.
113      * @param f Function to solve.
114      * @param min Lower bound for the interval.
115      * @param max Upper bound for the interval.
116      * @param allowedSolution The kind of solutions that the root-finding algorithm may
117      * accept as solutions.
118      * @return A value where the function is zero.
119      * @throws org.hipparchus.exception.MathIllegalArgumentException
120      * if the arguments do not satisfy the requirements specified by the solver.
121      * @throws org.hipparchus.exception.MathIllegalStateException if
122      * the allowed number of evaluations is exceeded.
123      */
124     T solve(int maxEval, CalculusFieldUnivariateFunction<T> f, T min, T max,
125             AllowedSolution allowedSolution);
126 
127     /**
128      * Solve for a zero in the given interval, start at {@code startValue}.
129      * A solver may require that the interval brackets a single zero root.
130      * Solvers that do require bracketing should be able to handle the case
131      * where one of the endpoints is itself a root.
132      *
133      * @param maxEval Maximum number of evaluations.
134      * @param f Function to solve.
135      * @param min Lower bound for the interval.
136      * @param max Upper bound for the interval.
137      * @param startValue Start value to use.
138      * @param allowedSolution The kind of solutions that the root-finding algorithm may
139      * accept as solutions.
140      * @return A value where the function is zero.
141      * @throws org.hipparchus.exception.MathIllegalArgumentException
142      * if the arguments do not satisfy the requirements specified by the solver.
143      * @throws org.hipparchus.exception.MathIllegalStateException if
144      * the allowed number of evaluations is exceeded.
145      */
146     T solve(int maxEval, CalculusFieldUnivariateFunction<T> f, T min, T max, T startValue,
147             AllowedSolution allowedSolution);
148 
149     /**
150      * Solve for a zero in the given interval and return a tolerance interval surrounding
151      * the root.
152      *
153      * <p> It is required that the starting interval brackets a root.
154      *
155      * @param maxEval Maximum number of evaluations.
156      * @param f       Function to solve.
157      * @param min     Lower bound for the interval. f(min) != 0.0.
158      * @param max     Upper bound for the interval. f(max) != 0.0.
159      * @return an interval [ta, tb] such that for some t in [ta, tb] f(t) == 0.0 or has a
160      * step wise discontinuity that crosses zero. Both end points also satisfy the
161      * convergence criteria so either one could be used as the root. That is the interval
162      * satisfies the condition (| tb - ta | &lt;= {@link #getAbsoluteAccuracy() absolute}
163      * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
164      * max(|f(ta)|, |f(tb)|) &lt;= {@link #getFunctionValueAccuracy()}) or there are no
165      * numbers in the field between ta and tb. The width of the interval (tb - ta) may be
166      * zero.
167      * @throws MathIllegalArgumentException if the arguments do not satisfy the
168      *                                      requirements specified by the solver.
169      * @throws MathIllegalStateException    if the allowed number of evaluations is
170      *                                      exceeded.
171      */
172     default Interval<T> solveInterval(int maxEval,
173                                       CalculusFieldUnivariateFunction<T> f,
174                                       T min,
175                                       T max)
176             throws MathIllegalArgumentException, MathIllegalStateException {
177         return this.solveInterval(maxEval, f, min, max, min.add(max.subtract(min).multiply(0.5)));
178     }
179 
180     /**
181      * Solve for a zero in the given interval and return a tolerance interval surrounding
182      * the root.
183      *
184      * <p> It is required that the starting interval brackets a root.
185      *
186      * @param maxEval    Maximum number of evaluations.
187      * @param startValue start value to use.
188      * @param f          Function to solve.
189      * @param min        Lower bound for the interval. f(min) != 0.0.
190      * @param max        Upper bound for the interval. f(max) != 0.0.
191      * @return an interval [ta, tb] such that for some t in [ta, tb] f(t) == 0.0 or has a
192      * step wise discontinuity that crosses zero. Both end points also satisfy the
193      * convergence criteria so either one could be used as the root. That is the interval
194      * satisfies the condition (| tb - ta | &lt;= {@link #getAbsoluteAccuracy() absolute}
195      * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
196      * max(|f(ta)|, |f(tb)|) &lt;= {@link #getFunctionValueAccuracy()}) or numbers in the
197      * field between ta and tb. The width of the interval (tb - ta) may be zero.
198      * @throws MathIllegalArgumentException if the arguments do not satisfy the
199      *                                      requirements specified by the solver.
200      * @throws MathIllegalStateException    if the allowed number of evaluations is
201      *                                      exceeded.
202      */
203     Interval<T> solveInterval(int maxEval,
204                               CalculusFieldUnivariateFunction<T> f,
205                               T min,
206                               T max,
207                               T startValue)
208             throws MathIllegalArgumentException, MathIllegalStateException;
209 
210     /**
211      * An interval of a function that brackets a root.
212      * <p>
213      * Contains two end points and the value of the function at the two end points.
214      *
215      * @see #solveInterval(int, CalculusFieldUnivariateFunction, CalculusFieldElement,
216      * CalculusFieldElement)
217      * @param <T> the element type
218      */
219     class Interval<T extends CalculusFieldElement<T>> {
220 
221         /** Abscissa on the left end of the interval. */
222         private final T leftAbscissa;
223         /** Function value at {@link #leftAbscissa}. */
224         private final T leftValue;
225         /** Abscissa on the right end of the interval, >= {@link #leftAbscissa}. */
226         private final T rightAbscissa;
227         /** Function value at {@link #rightAbscissa}. */
228         private final T rightValue;
229 
230         /**
231          * Construct a new interval with the given end points.
232          *
233          * @param leftAbscissa  is the abscissa value at the left side of the interval.
234          * @param leftValue     is the function value at {@code leftAbscissa}.
235          * @param rightAbscissa is the abscissa value on the right side of the interval.
236          *                      Must be greater than or equal to {@code leftAbscissa}.
237          * @param rightValue    is the function value at {@code rightAbscissa}.
238          */
239         public Interval(final T leftAbscissa,
240                         final T leftValue,
241                         final T rightAbscissa,
242                         final T rightValue) {
243             this.leftAbscissa = leftAbscissa;
244             this.leftValue = leftValue;
245             this.rightAbscissa = rightAbscissa;
246             this.rightValue = rightValue;
247         }
248 
249         /**
250          * Get the left abscissa.
251          *
252          * @return abscissa of the start of the interval.
253          */
254         public T getLeftAbscissa() {
255             return leftAbscissa;
256         }
257 
258         /**
259          * Get the right abscissa.
260          *
261          * @return abscissa of the end of the interval.
262          */
263         public T getRightAbscissa() {
264             return rightAbscissa;
265         }
266 
267         /**
268          * Get the function value at {@link #getLeftAbscissa()}.
269          *
270          * @return value of the function at the start of the interval.
271          */
272         public T getLeftValue() {
273             return leftValue;
274         }
275 
276         /**
277          * Get the function value at {@link #getRightAbscissa()}.
278          *
279          * @return value of the function at the end of the interval.
280          */
281         public T getRightValue() {
282             return rightValue;
283         }
284 
285         /**
286          * Get the abscissa corresponding to the allowed side.
287          *
288          * @param allowed side of the root.
289          * @return the abscissa on the selected side of the root.
290          */
291         public T getSide(final AllowedSolution allowed) {
292             final T xA = this.getLeftAbscissa();
293             final T yA = this.getLeftValue();
294             final T xB = this.getRightAbscissa();
295             switch (allowed) {
296                 case ANY_SIDE:
297                     final T absYA = this.getLeftValue().abs();
298                     final T absYB = this.getRightValue().abs();
299                     return absYA.subtract(absYB).getReal() < 0 ? xA : xB;
300                 case LEFT_SIDE:
301                     return xA;
302                 case RIGHT_SIDE:
303                     return xB;
304                 case BELOW_SIDE:
305                     return (yA.getReal() <= 0) ? xA : xB;
306                 case ABOVE_SIDE:
307                     return (yA.getReal() < 0) ? xB : xA;
308                 default:
309                     // this should never happen
310                     throw MathRuntimeException.createInternalError();
311             }
312         }
313 
314     }
315 
316 }