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 * http://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 package org.apache.commons.math.analysis.solvers; 18 19 import org.apache.commons.math.ConvergenceException; 20 import org.apache.commons.math.FunctionEvaluationException; 21 import org.apache.commons.math.MathRuntimeException; 22 import org.apache.commons.math.MaxIterationsExceededException; 23 import org.apache.commons.math.analysis.UnivariateRealFunction; 24 import org.apache.commons.math.exception.util.LocalizedFormats; 25 import org.apache.commons.math.util.FastMath; 26 27 28 /** 29 * Implements a modified version of the 30 * <a href="http://mathworld.wolfram.com/SecantMethod.html">secant method</a> 31 * for approximating a zero of a real univariate function. 32 * <p> 33 * The algorithm is modified to maintain bracketing of a root by successive 34 * approximations. Because of forced bracketing, convergence may be slower than 35 * the unrestricted secant algorithm. However, this implementation should in 36 * general outperform the 37 * <a href="http://mathworld.wolfram.com/MethodofFalsePosition.html"> 38 * regula falsi method.</a></p> 39 * <p> 40 * The function is assumed to be continuous but not necessarily smooth.</p> 41 * 42 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $ 43 */ 44 public class SecantSolver extends UnivariateRealSolverImpl { 45 46 /** 47 * Construct a solver for the given function. 48 * @param f function to solve. 49 * @deprecated as of 2.0 the function to solve is passed as an argument 50 * to the {@link #solve(UnivariateRealFunction, double, double)} or 51 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 52 * method. 53 */ 54 @Deprecated 55 public SecantSolver(UnivariateRealFunction f) { 56 super(f, 100, 1E-6); 57 } 58 59 /** 60 * Construct a solver. 61 * @deprecated in 2.2 (to be removed in 3.0). 62 */ 63 @Deprecated 64 public SecantSolver() { 65 super(100, 1E-6); 66 } 67 68 /** {@inheritDoc} */ 69 @Deprecated 70 public double solve(final double min, final double max) 71 throws ConvergenceException, FunctionEvaluationException { 72 return solve(f, min, max); 73 } 74 75 /** {@inheritDoc} */ 76 @Deprecated 77 public double solve(final double min, final double max, final double initial) 78 throws ConvergenceException, FunctionEvaluationException { 79 return solve(f, min, max, initial); 80 } 81 82 /** 83 * Find a zero in the given interval. 84 * 85 * @param f the function to solve 86 * @param min the lower bound for the interval 87 * @param max the upper bound for the interval 88 * @param initial the start value to use (ignored) 89 * @param maxEval Maximum number of evaluations. 90 * @return the value where the function is zero 91 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 92 * @throws FunctionEvaluationException if an error occurs evaluating the function 93 * @throws IllegalArgumentException if min is not less than max or the 94 * signs of the values of the function at the endpoints are not opposites 95 */ 96 @Override 97 public double solve(int maxEval, final UnivariateRealFunction f, 98 final double min, final double max, final double initial) 99 throws MaxIterationsExceededException, FunctionEvaluationException { 100 setMaximalIterationCount(maxEval); 101 return solve(f, min, max, initial); 102 } 103 104 /** 105 * Find a zero in the given interval. 106 * 107 * @param f the function to solve 108 * @param min the lower bound for the interval 109 * @param max the upper bound for the interval 110 * @param initial the start value to use (ignored) 111 * @return the value where the function is zero 112 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 113 * @throws FunctionEvaluationException if an error occurs evaluating the function 114 * @throws IllegalArgumentException if min is not less than max or the 115 * signs of the values of the function at the endpoints are not opposites 116 * @deprecated in 2.2 (to be removed in 3.0). 117 */ 118 @Deprecated 119 public double solve(final UnivariateRealFunction f, 120 final double min, final double max, final double initial) 121 throws MaxIterationsExceededException, FunctionEvaluationException { 122 return solve(f, min, max); 123 } 124 125 /** 126 * Find a zero in the given interval. 127 * @param f the function to solve 128 * @param min the lower bound for the interval. 129 * @param max the upper bound for the interval. 130 * @param maxEval Maximum number of evaluations. 131 * @return the value where the function is zero 132 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 133 * @throws FunctionEvaluationException if an error occurs evaluating the function 134 * @throws IllegalArgumentException if min is not less than max or the 135 * signs of the values of the function at the endpoints are not opposites 136 */ 137 @Override 138 public double solve(int maxEval, final UnivariateRealFunction f, 139 final double min, final double max) 140 throws MaxIterationsExceededException, FunctionEvaluationException { 141 setMaximalIterationCount(maxEval); 142 return solve(f, min, max); 143 } 144 145 /** 146 * Find a zero in the given interval. 147 * @param f the function to solve 148 * @param min the lower bound for the interval. 149 * @param max the upper bound for the interval. 150 * @return the value where the function is zero 151 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 152 * @throws FunctionEvaluationException if an error occurs evaluating the function 153 * @throws IllegalArgumentException if min is not less than max or the 154 * signs of the values of the function at the endpoints are not opposites 155 * @deprecated in 2.2 (to be removed in 3.0). 156 */ 157 @Deprecated 158 public double solve(final UnivariateRealFunction f, 159 final double min, final double max) 160 throws MaxIterationsExceededException, FunctionEvaluationException { 161 162 clearResult(); 163 verifyInterval(min, max); 164 165 // Index 0 is the old approximation for the root. 166 // Index 1 is the last calculated approximation for the root. 167 // Index 2 is a bracket for the root with respect to x0. 168 // OldDelta is the length of the bracketing interval of the last 169 // iteration. 170 double x0 = min; 171 double x1 = max; 172 double y0 = f.value(x0); 173 double y1 = f.value(x1); 174 175 // Verify bracketing 176 if (y0 * y1 >= 0) { 177 throw MathRuntimeException.createIllegalArgumentException( 178 LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, y0, y1); 179 } 180 181 double x2 = x0; 182 double y2 = y0; 183 double oldDelta = x2 - x1; 184 int i = 0; 185 while (i < maximalIterationCount) { 186 if (FastMath.abs(y2) < FastMath.abs(y1)) { 187 x0 = x1; 188 x1 = x2; 189 x2 = x0; 190 y0 = y1; 191 y1 = y2; 192 y2 = y0; 193 } 194 if (FastMath.abs(y1) <= functionValueAccuracy) { 195 setResult(x1, i); 196 return result; 197 } 198 if (FastMath.abs(oldDelta) < 199 FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy)) { 200 setResult(x1, i); 201 return result; 202 } 203 double delta; 204 if (FastMath.abs(y1) > FastMath.abs(y0)) { 205 // Function value increased in last iteration. Force bisection. 206 delta = 0.5 * oldDelta; 207 } else { 208 delta = (x0 - x1) / (1 - y0 / y1); 209 if (delta / oldDelta > 1) { 210 // New approximation falls outside bracket. 211 // Fall back to bisection. 212 delta = 0.5 * oldDelta; 213 } 214 } 215 x0 = x1; 216 y0 = y1; 217 x1 = x1 + delta; 218 y1 = f.value(x1); 219 if ((y1 > 0) == (y2 > 0)) { 220 // New bracket is (x0,x1). 221 x2 = x0; 222 y2 = y0; 223 } 224 oldDelta = x2 - x1; 225 i++; 226 } 227 throw new MaxIterationsExceededException(maximalIterationCount); 228 } 229 230 } 231
