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 }