001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math3.analysis.solvers; 019 020 import org.apache.commons.math3.util.FastMath; 021 import org.apache.commons.math3.analysis.UnivariateFunction; 022 import org.apache.commons.math3.exception.ConvergenceException; 023 import org.apache.commons.math3.exception.MathInternalError; 024 025 /** 026 * Base class for all bracketing <em>Secant</em>-based methods for root-finding 027 * (approximating a zero of a univariate real function). 028 * 029 * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and 030 * {@link IllinoisSolver <em>Illinois</em>} methods is based on the 031 * following article: M. Dowell and P. Jarratt, 032 * <em>A modified regula falsi method for computing the root of an 033 * equation</em>, BIT Numerical Mathematics, volume 11, number 2, 034 * pages 168-174, Springer, 1971.</p> 035 * 036 * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is 037 * based on the following article: M. Dowell and P. Jarratt, 038 * <em>The "Pegasus" method for computing the root of an equation</em>, 039 * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer, 040 * 1972.</p> 041 * 042 * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a 043 * bracketing method, so it is not implemented here. It has a separate 044 * implementation.</p> 045 * 046 * @since 3.0 047 * @version $Id: BaseSecantSolver.java 1379560 2012-08-31 19:40:30Z erans $ 048 */ 049 public abstract class BaseSecantSolver 050 extends AbstractUnivariateSolver 051 implements BracketedUnivariateSolver<UnivariateFunction> { 052 053 /** Default absolute accuracy. */ 054 protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; 055 056 /** The kinds of solutions that the algorithm may accept. */ 057 private AllowedSolution allowed; 058 059 /** The <em>Secant</em>-based root-finding method to use. */ 060 private final Method method; 061 062 /** 063 * Construct a solver. 064 * 065 * @param absoluteAccuracy Absolute accuracy. 066 * @param method <em>Secant</em>-based root-finding method to use. 067 */ 068 protected BaseSecantSolver(final double absoluteAccuracy, final Method method) { 069 super(absoluteAccuracy); 070 this.allowed = AllowedSolution.ANY_SIDE; 071 this.method = method; 072 } 073 074 /** 075 * Construct a solver. 076 * 077 * @param relativeAccuracy Relative accuracy. 078 * @param absoluteAccuracy Absolute accuracy. 079 * @param method <em>Secant</em>-based root-finding method to use. 080 */ 081 protected BaseSecantSolver(final double relativeAccuracy, 082 final double absoluteAccuracy, 083 final Method method) { 084 super(relativeAccuracy, absoluteAccuracy); 085 this.allowed = AllowedSolution.ANY_SIDE; 086 this.method = method; 087 } 088 089 /** 090 * Construct a solver. 091 * 092 * @param relativeAccuracy Maximum relative error. 093 * @param absoluteAccuracy Maximum absolute error. 094 * @param functionValueAccuracy Maximum function value error. 095 * @param method <em>Secant</em>-based root-finding method to use 096 */ 097 protected BaseSecantSolver(final double relativeAccuracy, 098 final double absoluteAccuracy, 099 final double functionValueAccuracy, 100 final Method method) { 101 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); 102 this.allowed = AllowedSolution.ANY_SIDE; 103 this.method = method; 104 } 105 106 /** {@inheritDoc} */ 107 public double solve(final int maxEval, final UnivariateFunction f, 108 final double min, final double max, 109 final AllowedSolution allowedSolution) { 110 return solve(maxEval, f, min, max, min + 0.5 * (max - min), allowedSolution); 111 } 112 113 /** {@inheritDoc} */ 114 public double solve(final int maxEval, final UnivariateFunction f, 115 final double min, final double max, final double startValue, 116 final AllowedSolution allowedSolution) { 117 this.allowed = allowedSolution; 118 return super.solve(maxEval, f, min, max, startValue); 119 } 120 121 /** {@inheritDoc} */ 122 @Override 123 public double solve(final int maxEval, final UnivariateFunction f, 124 final double min, final double max, final double startValue) { 125 return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE); 126 } 127 128 /** 129 * {@inheritDoc} 130 * 131 * @throws ConvergenceException if the algorithm failed due to finite 132 * precision. 133 */ 134 @Override 135 protected final double doSolve() 136 throws ConvergenceException, 137 MathInternalError { 138 // Get initial solution 139 double x0 = getMin(); 140 double x1 = getMax(); 141 double f0 = computeObjectiveValue(x0); 142 double f1 = computeObjectiveValue(x1); 143 144 // If one of the bounds is the exact root, return it. Since these are 145 // not under-approximations or over-approximations, we can return them 146 // regardless of the allowed solutions. 147 if (f0 == 0.0) { 148 return x0; 149 } 150 if (f1 == 0.0) { 151 return x1; 152 } 153 154 // Verify bracketing of initial solution. 155 verifyBracketing(x0, x1); 156 157 // Get accuracies. 158 final double ftol = getFunctionValueAccuracy(); 159 final double atol = getAbsoluteAccuracy(); 160 final double rtol = getRelativeAccuracy(); 161 162 // Keep track of inverted intervals, meaning that the left bound is 163 // larger than the right bound. 164 boolean inverted = false; 165 166 // Keep finding better approximations. 167 while (true) { 168 // Calculate the next approximation. 169 final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0)); 170 final double fx = computeObjectiveValue(x); 171 172 // If the new approximation is the exact root, return it. Since 173 // this is not an under-approximation or an over-approximation, 174 // we can return it regardless of the allowed solutions. 175 if (fx == 0.0) { 176 return x; 177 } 178 179 // Update the bounds with the new approximation. 180 if (f1 * fx < 0) { 181 // The value of x1 has switched to the other bound, thus inverting 182 // the interval. 183 x0 = x1; 184 f0 = f1; 185 inverted = !inverted; 186 } else { 187 switch (method) { 188 case ILLINOIS: 189 f0 *= 0.5; 190 break; 191 case PEGASUS: 192 f0 *= f1 / (f1 + fx); 193 break; 194 case REGULA_FALSI: 195 // Detect early that algorithm is stuck, instead of waiting 196 // for the maximum number of iterations to be exceeded. 197 if (x == x1) { 198 throw new ConvergenceException(); 199 } 200 break; 201 default: 202 // Should never happen. 203 throw new MathInternalError(); 204 } 205 } 206 // Update from [x0, x1] to [x0, x]. 207 x1 = x; 208 f1 = fx; 209 210 // If the function value of the last approximation is too small, 211 // given the function value accuracy, then we can't get closer to 212 // the root than we already are. 213 if (FastMath.abs(f1) <= ftol) { 214 switch (allowed) { 215 case ANY_SIDE: 216 return x1; 217 case LEFT_SIDE: 218 if (inverted) { 219 return x1; 220 } 221 break; 222 case RIGHT_SIDE: 223 if (!inverted) { 224 return x1; 225 } 226 break; 227 case BELOW_SIDE: 228 if (f1 <= 0) { 229 return x1; 230 } 231 break; 232 case ABOVE_SIDE: 233 if (f1 >= 0) { 234 return x1; 235 } 236 break; 237 default: 238 throw new MathInternalError(); 239 } 240 } 241 242 // If the current interval is within the given accuracies, we 243 // are satisfied with the current approximation. 244 if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), 245 atol)) { 246 switch (allowed) { 247 case ANY_SIDE: 248 return x1; 249 case LEFT_SIDE: 250 return inverted ? x1 : x0; 251 case RIGHT_SIDE: 252 return inverted ? x0 : x1; 253 case BELOW_SIDE: 254 return (f1 <= 0) ? x1 : x0; 255 case ABOVE_SIDE: 256 return (f1 >= 0) ? x1 : x0; 257 default: 258 throw new MathInternalError(); 259 } 260 } 261 } 262 } 263 264 /** <em>Secant</em>-based root-finding methods. */ 265 protected enum Method { 266 267 /** 268 * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or 269 * <em>False Position</em> method. 270 */ 271 REGULA_FALSI, 272 273 /** The {@link IllinoisSolver <em>Illinois</em>} method. */ 274 ILLINOIS, 275 276 /** The {@link PegasusSolver <em>Pegasus</em>} method. */ 277 PEGASUS; 278 279 } 280 }