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.analysis.UnivariateFunction;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.exception.MathIllegalStateException;
28  import org.hipparchus.exception.MathRuntimeException;
29  import org.hipparchus.util.FastMath;
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   * @param <F> Type of function to solve.
51   *
52   * @see AllowedSolution
53   */
54  public interface BracketedUnivariateSolver<F extends UnivariateFunction>
55      extends BaseUnivariateSolver<F> {
56  
57      /**
58       * Solve for a zero in the given interval.
59       * A solver may require that the interval brackets a single zero root.
60       * Solvers that do require bracketing should be able to handle the case
61       * where one of the endpoints is itself a root.
62       *
63       * @param maxEval Maximum number of evaluations.
64       * @param f Function to solve.
65       * @param min Lower bound for the interval.
66       * @param max Upper bound for the interval.
67       * @param allowedSolution The kind of solutions that the root-finding algorithm may
68       * accept as solutions.
69       * @return A value where the function is zero.
70       * @throws org.hipparchus.exception.MathIllegalArgumentException
71       * if the arguments do not satisfy the requirements specified by the solver.
72       * @throws org.hipparchus.exception.MathIllegalStateException if
73       * the allowed number of evaluations is exceeded.
74       */
75      double solve(int maxEval, F f, double min, double max,
76                   AllowedSolution allowedSolution);
77  
78      /**
79       * Solve for a zero in the given interval, start at {@code startValue}.
80       * A solver may require that the interval brackets a single zero root.
81       * Solvers that do require bracketing should be able to handle the case
82       * where one of the endpoints is itself a root.
83       *
84       * @param maxEval Maximum number of evaluations.
85       * @param f Function to solve.
86       * @param min Lower bound for the interval.
87       * @param max Upper bound for the interval.
88       * @param startValue Start value to use.
89       * @param allowedSolution The kind of solutions that the root-finding algorithm may
90       * accept as solutions.
91       * @return A value where the function is zero.
92       * @throws org.hipparchus.exception.MathIllegalArgumentException
93       * if the arguments do not satisfy the requirements specified by the solver.
94       * @throws org.hipparchus.exception.MathIllegalStateException if
95       * the allowed number of evaluations is exceeded.
96       */
97      double solve(int maxEval, F f, double min, double max, double startValue,
98                   AllowedSolution allowedSolution);
99  
100     /**
101      * Solve for a zero in the given interval and return a tolerance interval surrounding
102      * the root.
103      *
104      * <p> It is required that the starting interval brackets a root or that the function
105      * value at either end point is 0.0.
106      *
107      * @param maxEval Maximum number of evaluations.
108      * @param f       Function to solve.
109      * @param min     Lower bound for the interval.
110      * @param max     Upper bound for the interval. Must be greater than {@code min}.
111      * @return an interval [ta, tb] such that for some t in [ta, tb] f(t) == 0.0 or has a
112      * step wise discontinuity that crosses zero. Both end points also satisfy the
113      * convergence criteria so either one could be used as the root. That is the interval
114      * satisfies the condition (| tb - ta | &lt;= {@link #getAbsoluteAccuracy() absolute}
115      * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
116      * max(|f(ta)|, |f(tb)|) &lt;= {@link #getFunctionValueAccuracy()}) or there are no
117      * floating point numbers between ta and tb. The width of the interval (tb - ta) may
118      * be zero.
119      * @throws MathIllegalArgumentException if the arguments do not satisfy the
120      *                                      requirements specified by the solver.
121      * @throws MathIllegalStateException    if the allowed number of evaluations is
122      *                                      exceeded.
123      */
124     default Interval solveInterval(int maxEval, F f, double min, double max)
125             throws MathIllegalArgumentException, MathIllegalStateException {
126         return this.solveInterval(maxEval, f, min, max, min + 0.5 * (max - min));
127     }
128 
129     /**
130      * Solve for a zero in the given interval and return a tolerance interval surrounding
131      * the root.
132      *
133      * <p> It is required that the starting interval brackets a root or that the function
134      * value at either end point is 0.0.
135      *
136      * @param maxEval    Maximum number of evaluations.
137      * @param startValue start value to use. Must be in the interval [min, max].
138      * @param f          Function to solve.
139      * @param min        Lower bound for the interval.
140      * @param max     Upper bound for the interval. Must be greater than {@code min}.
141      * @return an interval [ta, tb] such that for some t in [ta, tb] f(t) == 0.0 or has a
142      * step wise discontinuity that crosses zero. Both end points also satisfy the
143      * convergence criteria so either one could be used as the root. That is the interval
144      * satisfies the condition (| tb - ta | &lt;= {@link #getAbsoluteAccuracy() absolute}
145      * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
146      * max(|f(ta)|, |f(tb)|) &lt;= {@link #getFunctionValueAccuracy()}) or there are no
147      * floating point numbers between ta and tb. The width of the interval (tb - ta) may
148      * be zero.
149      * @throws MathIllegalArgumentException if the arguments do not satisfy the
150      *                                      requirements specified by the solver.
151      * @throws MathIllegalStateException    if the allowed number of evaluations is
152      *                                      exceeded.
153      */
154     Interval solveInterval(int maxEval, F f, double min, double max, double startValue)
155             throws MathIllegalArgumentException, MathIllegalStateException;
156 
157     /**
158      * An interval of a function that brackets a root.
159      *
160      * <p> Contains two end points and the value of the function at the two end points.
161      *
162      * @see #solveInterval(int, UnivariateFunction, double, double, double)
163      */
164     class Interval {
165 
166         /** Abscissa on the left end of the interval. */
167         private final double leftAbscissa;
168         /** Function value at {@link #leftAbscissa}. */
169         private final double leftValue;
170         /** Abscissa on the right end of the interval, >= {@link #leftAbscissa}. */
171         private final double rightAbscissa;
172         /** Function value at {@link #rightAbscissa}. */
173         private final double rightValue;
174 
175         /**
176          * Construct a new interval with the given end points.
177          *
178          * @param leftAbscissa  is the abscissa value at the left side of the interval.
179          * @param leftValue     is the function value at {@code leftAbscissa}.
180          * @param rightAbscissa is the abscissa value on the right side of the interval.
181          *                      Must be greater than or equal to {@code leftAbscissa}.
182          * @param rightValue    is the function value at {@code rightAbscissa}.
183          */
184         public Interval(final double leftAbscissa,
185                         final double leftValue,
186                         final double rightAbscissa,
187                         final double rightValue) {
188             this.leftAbscissa = leftAbscissa;
189             this.leftValue = leftValue;
190             this.rightAbscissa = rightAbscissa;
191             this.rightValue = rightValue;
192         }
193 
194         /**
195          * Get the left abscissa.
196          *
197          * @return abscissa of the start of the interval.
198          */
199         public double getLeftAbscissa() {
200             return leftAbscissa;
201         }
202 
203         /**
204          * Get the right abscissa.
205          *
206          * @return abscissa of the end of the interval.
207          */
208         public double getRightAbscissa() {
209             return rightAbscissa;
210         }
211 
212         /**
213          * Get the function value at {@link #getLeftAbscissa()}.
214          *
215          * @return value of the function at the start of the interval.
216          */
217         public double getLeftValue() {
218             return leftValue;
219         }
220 
221         /**
222          * Get the function value at {@link #getRightAbscissa()}.
223          *
224          * @return value of the function at the end of the interval.
225          */
226         public double getRightValue() {
227             return rightValue;
228         }
229 
230         /**
231          * Get the abscissa corresponding to the allowed side.
232          *
233          * @param allowed side of the root.
234          * @return the abscissa on the selected side of the root.
235          */
236         public double getSide(final AllowedSolution allowed) {
237             final double xA = this.getLeftAbscissa();
238             final double yA = this.getLeftValue();
239             final double xB = this.getRightAbscissa();
240             switch (allowed) {
241                 case ANY_SIDE:
242                     final double absYA = FastMath.abs(this.getLeftValue());
243                     final double absYB = FastMath.abs(this.getRightValue());
244                     return absYA < absYB ? xA : xB;
245                 case LEFT_SIDE:
246                     return xA;
247                 case RIGHT_SIDE:
248                     return xB;
249                 case BELOW_SIDE:
250                     return (yA <= 0) ? xA : xB;
251                 case ABOVE_SIDE:
252                     return (yA < 0) ? xB : xA;
253                 default:
254                     // this should never happen
255                     throw MathRuntimeException.createInternalError();
256             }
257         }
258 
259     }
260 
261 }