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 package org.apache.commons.math.analysis.solvers; 018 019 import org.apache.commons.math.ConvergenceException; 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.MathRuntimeException; 022 import org.apache.commons.math.MaxIterationsExceededException; 023 import org.apache.commons.math.analysis.UnivariateRealFunction; 024 import org.apache.commons.math.exception.util.LocalizedFormats; 025 import org.apache.commons.math.util.FastMath; 026 027 028 /** 029 * Implements a modified version of the 030 * <a href="http://mathworld.wolfram.com/SecantMethod.html">secant method</a> 031 * for approximating a zero of a real univariate function. 032 * <p> 033 * The algorithm is modified to maintain bracketing of a root by successive 034 * approximations. Because of forced bracketing, convergence may be slower than 035 * the unrestricted secant algorithm. However, this implementation should in 036 * general outperform the 037 * <a href="http://mathworld.wolfram.com/MethodofFalsePosition.html"> 038 * regula falsi method.</a></p> 039 * <p> 040 * The function is assumed to be continuous but not necessarily smooth.</p> 041 * 042 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $ 043 */ 044 public class SecantSolver extends UnivariateRealSolverImpl { 045 046 /** 047 * Construct a solver for the given function. 048 * @param f function to solve. 049 * @deprecated as of 2.0 the function to solve is passed as an argument 050 * to the {@link #solve(UnivariateRealFunction, double, double)} or 051 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 052 * method. 053 */ 054 @Deprecated 055 public SecantSolver(UnivariateRealFunction f) { 056 super(f, 100, 1E-6); 057 } 058 059 /** 060 * Construct a solver. 061 * @deprecated in 2.2 (to be removed in 3.0). 062 */ 063 @Deprecated 064 public SecantSolver() { 065 super(100, 1E-6); 066 } 067 068 /** {@inheritDoc} */ 069 @Deprecated 070 public double solve(final double min, final double max) 071 throws ConvergenceException, FunctionEvaluationException { 072 return solve(f, min, max); 073 } 074 075 /** {@inheritDoc} */ 076 @Deprecated 077 public double solve(final double min, final double max, final double initial) 078 throws ConvergenceException, FunctionEvaluationException { 079 return solve(f, min, max, initial); 080 } 081 082 /** 083 * Find a zero in the given interval. 084 * 085 * @param f the function to solve 086 * @param min the lower bound for the interval 087 * @param max the upper bound for the interval 088 * @param initial the start value to use (ignored) 089 * @param maxEval Maximum number of evaluations. 090 * @return the value where the function is zero 091 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 092 * @throws FunctionEvaluationException if an error occurs evaluating the function 093 * @throws IllegalArgumentException if min is not less than max or the 094 * signs of the values of the function at the endpoints are not opposites 095 */ 096 @Override 097 public double solve(int maxEval, final UnivariateRealFunction f, 098 final double min, final double max, final double initial) 099 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 }