1 /* 2 * Copyright 2003-2004 The Apache Software Foundation. 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 package org.apache.commons.math.analysis; 17 18 import java.io.Serializable; 19 20 import org.apache.commons.math.ConvergenceException; 21 import org.apache.commons.math.FunctionEvaluationException; 22 23 24 /** 25 * Implements a modified version of the 26 * <a href="http://mathworld.wolfram.com/SecantMethod.html">secant method</a> 27 * for approximating a zero of a real univariate function. 28 * <p> 29 * The algorithm is modified to maintain bracketing of a root by successive 30 * approximations. Because of forced bracketing, convergence may be slower than 31 * the unrestricted secant algorithm. However, this implementation should in 32 * general outperform the 33 * <a href="http://mathworld.wolfram.com/MethodofFalsePosition.html"> 34 * regula falsi method.</a> 35 * <p> 36 * The function is assumed to be continuous but not necessarily smooth. 37 * 38 * @version $Revision: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $ 39 */ 40 public class SecantSolver extends UnivariateRealSolverImpl implements Serializable { 41 42 /** Serializable version identifier */ 43 private static final long serialVersionUID = 1984971194738974867L; 44 45 /** 46 * Construct a solver for the given function. 47 * @param f function to solve. 48 */ 49 public SecantSolver(UnivariateRealFunction f) { 50 super(f, 100, 1E-6); 51 } 52 53 /** 54 * Find a zero in the given interval. 55 * 56 * @param min the lower bound for the interval 57 * @param max the upper bound for the interval 58 * @param initial the start value to use (ignored) 59 * @return the value where the function is zero 60 * @throws ConvergenceException if the maximum iteration count is exceeded 61 * @throws FunctionEvaluationException if an error occurs evaluating the 62 * function 63 * @throws IllegalArgumentException if min is not less than max or the 64 * signs of the values of the function at the endpoints are not opposites 65 */ 66 public double solve(double min, double max, double initial) 67 throws ConvergenceException, FunctionEvaluationException { 68 69 return solve(min, max); 70 } 71 72 /** 73 * Find a zero in the given interval. 74 * @param min the lower bound for the interval. 75 * @param max the upper bound for the interval. 76 * @return the value where the function is zero 77 * @throws ConvergenceException if the maximum iteration count is exceeded 78 * @throws FunctionEvaluationException if an error occurs evaluating the 79 * function 80 * @throws IllegalArgumentException if min is not less than max or the 81 * signs of the values of the function at the endpoints are not opposites 82 */ 83 public double solve(double min, double max) throws ConvergenceException, 84 FunctionEvaluationException { 85 86 clearResult(); 87 verifyInterval(min, max); 88 89 // Index 0 is the old approximation for the root. 90 // Index 1 is the last calculated approximation for the root. 91 // Index 2 is a bracket for the root with respect to x0. 92 // OldDelta is the length of the bracketing interval of the last 93 // iteration. 94 double x0 = min; 95 double x1 = max; 96 double y0 = f.value(x0); 97 double y1 = f.value(x1); 98 99 // Verify bracketing 100 if (y0 * y1 >= 0) { 101 throw new IllegalArgumentException 102 ("Function values at endpoints do not have different signs." + 103 " Endpoints: [" + min + "," + max + "]" + 104 " Values: [" + y0 + "," + y1 + "]"); 105 } 106 107 double x2 = x0; 108 double y2 = y0; 109 double oldDelta = x2 - x1; 110 int i = 0; 111 while (i < maximalIterationCount) { 112 if (Math.abs(y2) < Math.abs(y1)) { 113 x0 = x1; 114 x1 = x2; 115 x2 = x0; 116 y0 = y1; 117 y1 = y2; 118 y2 = y0; 119 } 120 if (Math.abs(y1) <= functionValueAccuracy) { 121 setResult(x1, i); 122 return result; 123 } 124 if (Math.abs(oldDelta) < 125 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy)) { 126 setResult(x1, i); 127 return result; 128 } 129 double delta; 130 if (Math.abs(y1) > Math.abs(y0)) { 131 // Function value increased in last iteration. Force bisection. 132 delta = 0.5 * oldDelta; 133 } else { 134 delta = (x0 - x1) / (1 - y0 / y1); 135 if (delta / oldDelta > 1) { 136 // New approximation falls outside bracket. 137 // Fall back to bisection. 138 delta = 0.5 * oldDelta; 139 } 140 } 141 x0 = x1; 142 y0 = y1; 143 x1 = x1 + delta; 144 y1 = f.value(x1); 145 if ((y1 > 0) == (y2 > 0)) { 146 // New bracket is (x0,x1). 147 x2 = x0; 148 y2 = y0; 149 } 150 oldDelta = x2 - x1; 151 i++; 152 } 153 throw new ConvergenceException("Maximal iteration number exceeded" + i); 154 } 155 156 }