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.math3.analysis.solvers; 018 019 020 import org.apache.commons.math3.exception.NoBracketingException; 021 import org.apache.commons.math3.exception.TooManyEvaluationsException; 022 import org.apache.commons.math3.exception.NumberIsTooLargeException; 023 import org.apache.commons.math3.util.FastMath; 024 import org.apache.commons.math3.util.Precision; 025 026 /** 027 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html"> 028 * Brent algorithm</a> for finding zeros of real univariate functions. 029 * The function should be continuous but not necessarily smooth. 030 * The {@code solve} method returns a zero {@code x} of the function {@code f} 031 * in the given interval {@code [a, b]} to within a tolerance 032 * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and 033 * {@code t} is the absolute accuracy. 034 * The given interval must bracket the root. 035 * 036 * @version $Id: BrentSolver.java 1379560 2012-08-31 19:40:30Z erans $ 037 */ 038 public class BrentSolver extends AbstractUnivariateSolver { 039 040 /** Default absolute accuracy. */ 041 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; 042 043 /** 044 * Construct a solver with default accuracy (1e-6). 045 */ 046 public BrentSolver() { 047 this(DEFAULT_ABSOLUTE_ACCURACY); 048 } 049 /** 050 * Construct a solver. 051 * 052 * @param absoluteAccuracy Absolute accuracy. 053 */ 054 public BrentSolver(double absoluteAccuracy) { 055 super(absoluteAccuracy); 056 } 057 /** 058 * Construct a solver. 059 * 060 * @param relativeAccuracy Relative accuracy. 061 * @param absoluteAccuracy Absolute accuracy. 062 */ 063 public BrentSolver(double relativeAccuracy, 064 double absoluteAccuracy) { 065 super(relativeAccuracy, absoluteAccuracy); 066 } 067 /** 068 * Construct a solver. 069 * 070 * @param relativeAccuracy Relative accuracy. 071 * @param absoluteAccuracy Absolute accuracy. 072 * @param functionValueAccuracy Function value accuracy. 073 */ 074 public BrentSolver(double relativeAccuracy, 075 double absoluteAccuracy, 076 double functionValueAccuracy) { 077 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); 078 } 079 080 /** 081 * {@inheritDoc} 082 */ 083 @Override 084 protected double doSolve() 085 throws NoBracketingException, 086 TooManyEvaluationsException, 087 NumberIsTooLargeException { 088 double min = getMin(); 089 double max = getMax(); 090 final double initial = getStartValue(); 091 final double functionValueAccuracy = getFunctionValueAccuracy(); 092 093 verifySequence(min, initial, max); 094 095 // Return the initial guess if it is good enough. 096 double yInitial = computeObjectiveValue(initial); 097 if (FastMath.abs(yInitial) <= functionValueAccuracy) { 098 return initial; 099 } 100 101 // Return the first endpoint if it is good enough. 102 double yMin = computeObjectiveValue(min); 103 if (FastMath.abs(yMin) <= functionValueAccuracy) { 104 return min; 105 } 106 107 // Reduce interval if min and initial bracket the root. 108 if (yInitial * yMin < 0) { 109 return brent(min, initial, yMin, yInitial); 110 } 111 112 // Return the second endpoint if it is good enough. 113 double yMax = computeObjectiveValue(max); 114 if (FastMath.abs(yMax) <= functionValueAccuracy) { 115 return max; 116 } 117 118 // Reduce interval if initial and max bracket the root. 119 if (yInitial * yMax < 0) { 120 return brent(initial, max, yInitial, yMax); 121 } 122 123 throw new NoBracketingException(min, max, yMin, yMax); 124 } 125 126 /** 127 * Search for a zero inside the provided interval. 128 * This implementation is based on the algorithm described at page 58 of 129 * the book 130 * <quote> 131 * <b>Algorithms for Minimization Without Derivatives</b> 132 * <it>Richard P. Brent</it> 133 * Dover 0-486-41998-3 134 * </quote> 135 * 136 * @param lo Lower bound of the search interval. 137 * @param hi Higher bound of the search interval. 138 * @param fLo Function value at the lower bound of the search interval. 139 * @param fHi Function value at the higher bound of the search interval. 140 * @return the value where the function is zero. 141 */ 142 private double brent(double lo, double hi, 143 double fLo, double fHi) { 144 double a = lo; 145 double fa = fLo; 146 double b = hi; 147 double fb = fHi; 148 double c = a; 149 double fc = fa; 150 double d = b - a; 151 double e = d; 152 153 final double t = getAbsoluteAccuracy(); 154 final double eps = getRelativeAccuracy(); 155 156 while (true) { 157 if (FastMath.abs(fc) < FastMath.abs(fb)) { 158 a = b; 159 b = c; 160 c = a; 161 fa = fb; 162 fb = fc; 163 fc = fa; 164 } 165 166 final double tol = 2 * eps * FastMath.abs(b) + t; 167 final double m = 0.5 * (c - b); 168 169 if (FastMath.abs(m) <= tol || 170 Precision.equals(fb, 0)) { 171 return b; 172 } 173 if (FastMath.abs(e) < tol || 174 FastMath.abs(fa) <= FastMath.abs(fb)) { 175 // Force bisection. 176 d = m; 177 e = d; 178 } else { 179 double s = fb / fa; 180 double p; 181 double q; 182 // The equality test (a == c) is intentional, 183 // it is part of the original Brent's method and 184 // it should NOT be replaced by proximity test. 185 if (a == c) { 186 // Linear interpolation. 187 p = 2 * m * s; 188 q = 1 - s; 189 } else { 190 // Inverse quadratic interpolation. 191 q = fa / fc; 192 final double r = fb / fc; 193 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); 194 q = (q - 1) * (r - 1) * (s - 1); 195 } 196 if (p > 0) { 197 q = -q; 198 } else { 199 p = -p; 200 } 201 s = e; 202 e = d; 203 if (p >= 1.5 * m * q - FastMath.abs(tol * q) || 204 p >= FastMath.abs(0.5 * s * q)) { 205 // Inverse quadratic interpolation gives a value 206 // in the wrong direction, or progress is slow. 207 // Fall back to bisection. 208 d = m; 209 e = d; 210 } else { 211 d = p / q; 212 } 213 } 214 a = b; 215 fa = fb; 216 217 if (FastMath.abs(d) > tol) { 218 b += d; 219 } else if (m > 0) { 220 b += tol; 221 } else { 222 b -= tol; 223 } 224 fb = computeObjectiveValue(b); 225 if ((fb > 0 && fc > 0) || 226 (fb <= 0 && fc <= 0)) { 227 c = a; 228 fc = fa; 229 d = b - a; 230 e = d; 231 } 232 } 233 } 234 }