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 ε 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} - ε, {@code v} + ε).
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 ρ 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} - ρ {@code v}, {@code v} + ρ {@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 | <= {@link #getAbsoluteAccuracy() absolute}
163 * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
164 * max(|f(ta)|, |f(tb)|) <= {@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 | <= {@link #getAbsoluteAccuracy() absolute}
195 * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or (
196 * max(|f(ta)|, |f(tb)|) <= {@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 }