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