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.numbers.rootfinder; 018 019import java.util.function.DoubleUnaryOperator; 020 021/** 022 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html"> 023 * Brent algorithm</a> for finding zeros of real univariate functions. 024 * The function should be continuous but not necessarily smooth. 025 * The {@code solve} method returns a zero {@code x} of the function {@code f} 026 * in the given interval {@code [a, b]} to within a tolerance 027 * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and 028 * {@code t} is the absolute accuracy. 029 * <p>The given interval must bracket the root.</p> 030 * <p> 031 * The reference implementation is given in chapter 4 of 032 * <blockquote> 033 * <b>Algorithms for Minimization Without Derivatives</b>, 034 * <em>Richard P. Brent</em>, 035 * Dover, 2002 036 * </blockquote> 037 */ 038public class BrentSolver { 039 /** Relative accuracy. */ 040 private final double relativeAccuracy; 041 /** Absolute accuracy. */ 042 private final double absoluteAccuracy; 043 /** Function accuracy. */ 044 private final double functionValueAccuracy; 045 046 /** 047 * Construct a solver. 048 * 049 * @param relativeAccuracy Relative accuracy. 050 * @param absoluteAccuracy Absolute accuracy. 051 * @param functionValueAccuracy Function value accuracy. 052 */ 053 public BrentSolver(double relativeAccuracy, 054 double absoluteAccuracy, 055 double functionValueAccuracy) { 056 this.relativeAccuracy = relativeAccuracy; 057 this.absoluteAccuracy = absoluteAccuracy; 058 this.functionValueAccuracy = functionValueAccuracy; 059 } 060 061 /** 062 * Search the function's zero within the given interval. 063 * 064 * @param func Function to solve. 065 * @param min Lower bound. 066 * @param max Upper bound. 067 * @return the root. 068 * @throws IllegalArgumentException if {@code min > max}. 069 * @throws IllegalArgumentException if the given interval does 070 * not bracket the root. 071 */ 072 public double findRoot(DoubleUnaryOperator func, 073 double min, 074 double max) { 075 // Avoid overflow computing the initial value: 0.5 * (min + max) 076 // Note: This sum is invalid if min == max == Double.MIN_VALUE 077 // so detect this edge case. It will raise a bracketing exception 078 // if min is not the root within the configured function accuracy; 079 // otherwise min is returned. 080 final double initial = min == max ? min : 0.5 * min + 0.5 * max; 081 return findRoot(func, min, initial, max); 082 } 083 084 /** 085 * Search the function's zero within the given interval, 086 * starting from the given estimate. 087 * 088 * @param func Function to solve. 089 * @param min Lower bound. 090 * @param initial Initial guess. 091 * @param max Upper bound. 092 * @return the root. 093 * @throws IllegalArgumentException if {@code min > max} or 094 * {@code initial} is not in the {@code [min, max]} interval. 095 * @throws IllegalArgumentException if the given interval does 096 * not bracket the root. 097 */ 098 public double findRoot(DoubleUnaryOperator func, 099 double min, 100 double initial, 101 double max) { 102 if (min > max) { 103 throw new SolverException(SolverException.TOO_LARGE, min, max); 104 } 105 if (initial < min || 106 initial > max) { 107 throw new SolverException(SolverException.OUT_OF_RANGE, initial, min, max); 108 } 109 110 // Return the initial guess if it is good enough. 111 final double yInitial = func.applyAsDouble(initial); 112 if (Math.abs(yInitial) <= functionValueAccuracy) { 113 return initial; 114 } 115 116 // Return the first endpoint if it is good enough. 117 final double yMin = func.applyAsDouble(min); 118 if (Math.abs(yMin) <= functionValueAccuracy) { 119 return min; 120 } 121 122 // Reduce interval if min and initial bracket the root. 123 if (Double.compare(yInitial * yMin, 0.0) < 0) { 124 return brent(func, min, initial, yMin, yInitial); 125 } 126 127 // Return the second endpoint if it is good enough. 128 final double yMax = func.applyAsDouble(max); 129 if (Math.abs(yMax) <= functionValueAccuracy) { 130 return max; 131 } 132 133 // Reduce interval if initial and max bracket the root. 134 if (Double.compare(yInitial * yMax, 0.0) < 0) { 135 return brent(func, initial, max, yInitial, yMax); 136 } 137 138 throw new SolverException(SolverException.BRACKETING, min, yMin, max, yMax); 139 } 140 141 /** 142 * Search for a zero inside the provided interval. 143 * This implementation is based on the algorithm described at page 58 of 144 * the book 145 * <blockquote> 146 * <b>Algorithms for Minimization Without Derivatives</b>, 147 * <i>Richard P. Brent</i>, 148 * Dover 0-486-41998-3 149 * </blockquote> 150 * 151 * @param func Function to solve. 152 * @param lo Lower bound of the search interval. 153 * @param hi Higher bound of the search interval. 154 * @param fLo Function value at the lower bound of the search interval. 155 * @param fHi Function value at the higher bound of the search interval. 156 * @return the value where the function is zero. 157 */ 158 private double brent(DoubleUnaryOperator func, 159 double lo, double hi, 160 double fLo, double fHi) { 161 double a = lo; 162 double fa = fLo; 163 double b = hi; 164 double fb = fHi; 165 double c = a; 166 double fc = fa; 167 double d = b - a; 168 double e = d; 169 170 final double t = absoluteAccuracy; 171 final double eps = relativeAccuracy; 172 173 while (true) { 174 if (Math.abs(fc) < Math.abs(fb)) { 175 a = b; 176 b = c; 177 c = a; 178 fa = fb; 179 fb = fc; 180 fc = fa; 181 } 182 183 final double tol = 2 * eps * Math.abs(b) + t; 184 final double m = 0.5 * (c - b); 185 186 if (Math.abs(m) <= tol || 187 equalsZero(fb)) { 188 return b; 189 } 190 if (Math.abs(e) < tol || 191 Math.abs(fa) <= Math.abs(fb)) { 192 // Force bisection. 193 d = m; 194 e = d; 195 } else { 196 final double s = fb / fa; 197 double p; 198 double q; 199 // The equality test (a == c) is intentional, 200 // it is part of the original Brent's method and 201 // it should NOT be replaced by proximity test. 202 if (a == c) { 203 // Linear interpolation. 204 p = 2 * m * s; 205 q = 1 - s; 206 } else { 207 // Inverse quadratic interpolation. 208 q = fa / fc; 209 final double r = fb / fc; 210 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); 211 q = (q - 1) * (r - 1) * (s - 1); 212 } 213 if (p > 0) { 214 q = -q; 215 } else { 216 p = -p; 217 } 218 if (p >= 1.5 * m * q - Math.abs(tol * q) || 219 p >= Math.abs(0.5 * e * q)) { 220 // Inverse quadratic interpolation gives a value 221 // in the wrong direction, or progress is slow. 222 // Fall back to bisection. 223 d = m; 224 e = d; 225 } else { 226 e = d; 227 d = p / q; 228 } 229 } 230 a = b; 231 fa = fb; 232 233 if (Math.abs(d) > tol) { 234 b += d; 235 } else if (m > 0) { 236 b += tol; 237 } else { 238 b -= tol; 239 } 240 fb = func.applyAsDouble(b); 241 if ((fb > 0 && fc > 0) || 242 (fb <= 0 && fc <= 0)) { 243 c = a; 244 fc = fa; 245 d = b - a; 246 e = d; 247 } 248 } 249 } 250 251 /** 252 * Return true if the value is within 1 ULP of zero. 253 * 254 * @param value Value 255 * @return true if zero within a 1 ULP tolerance 256 */ 257 private static boolean equalsZero(double value) { 258 return Math.abs(value) <= Double.MIN_VALUE; 259 } 260}