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