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 | <= {@link #getAbsoluteAccuracy() absolute} 115 * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or ( 116 * max(|f(ta)|, |f(tb)|) <= {@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 | <= {@link #getAbsoluteAccuracy() absolute} 145 * accuracy + max(ta, tb) * {@link #getRelativeAccuracy() relative} accuracy) or ( 146 * max(|f(ta)|, |f(tb)|) <= {@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 }