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