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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.solvers;
23
24
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.exception.MathIllegalStateException;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.Precision;
30
31 /**
32 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
33 * Brent algorithm</a> for finding zeros of real univariate functions.
34 * The function should be continuous but not necessarily smooth.
35 * The {@code solve} method returns a zero {@code x} of the function {@code f}
36 * in the given interval {@code [a, b]} to within a tolerance
37 * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
38 * {@code t} is the absolute accuracy.
39 * <p>The given interval must bracket the root.</p>
40 * <p>
41 * The reference implementation is given in chapter 4 of
42 * <blockquote>
43 * <b>Algorithms for Minimization Without Derivatives</b>,
44 * <em>Richard P. Brent</em>,
45 * Dover, 2002
46 * </blockquote>
47 *
48 * @see BaseAbstractUnivariateSolver
49 */
50 public class BrentSolver extends AbstractUnivariateSolver {
51
52 /** Default absolute accuracy. */
53 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
54
55 /**
56 * Construct a solver with default absolute accuracy (1e-6).
57 */
58 public BrentSolver() {
59 this(DEFAULT_ABSOLUTE_ACCURACY);
60 }
61 /**
62 * Construct a solver.
63 *
64 * @param absoluteAccuracy Absolute accuracy.
65 */
66 public BrentSolver(double absoluteAccuracy) {
67 super(absoluteAccuracy);
68 }
69 /**
70 * Construct a solver.
71 *
72 * @param relativeAccuracy Relative accuracy.
73 * @param absoluteAccuracy Absolute accuracy.
74 */
75 public BrentSolver(double relativeAccuracy,
76 double absoluteAccuracy) {
77 super(relativeAccuracy, absoluteAccuracy);
78 }
79 /**
80 * Construct a solver.
81 *
82 * @param relativeAccuracy Relative accuracy.
83 * @param absoluteAccuracy Absolute accuracy.
84 * @param functionValueAccuracy Function value accuracy.
85 *
86 * @see BaseAbstractUnivariateSolver#BaseAbstractUnivariateSolver(double,double,double)
87 */
88 public BrentSolver(double relativeAccuracy,
89 double absoluteAccuracy,
90 double functionValueAccuracy) {
91 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
92 }
93
94 /**
95 * {@inheritDoc}
96 */
97 @Override
98 protected double doSolve()
99 throws MathIllegalArgumentException, MathIllegalStateException {
100 double min = getMin();
101 double max = getMax();
102 final double initial = getStartValue();
103 final double functionValueAccuracy = getFunctionValueAccuracy();
104
105 verifySequence(min, initial, max);
106
107 // Return the initial guess if it is good enough.
108 double yInitial = computeObjectiveValue(initial);
109 if (FastMath.abs(yInitial) <= functionValueAccuracy) {
110 return initial;
111 }
112
113 // Return the first endpoint if it is good enough.
114 double yMin = computeObjectiveValue(min);
115 if (FastMath.abs(yMin) <= functionValueAccuracy) {
116 return min;
117 }
118
119 // Reduce interval if min and initial bracket the root.
120 if (yInitial * yMin < 0) {
121 return brent(min, initial, yMin, yInitial);
122 }
123
124 // Return the second endpoint if it is good enough.
125 double yMax = computeObjectiveValue(max);
126 if (FastMath.abs(yMax) <= functionValueAccuracy) {
127 return max;
128 }
129
130 // Reduce interval if initial and max bracket the root.
131 if (yInitial * yMax < 0) {
132 return brent(initial, max, yInitial, yMax);
133 }
134
135 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
136 min, max, yMin, yMax);
137 }
138
139 /**
140 * Search for a zero inside the provided interval.
141 * This implementation is based on the algorithm described at page 58 of
142 * the book
143 * <blockquote>
144 * <b>Algorithms for Minimization Without Derivatives</b>,
145 * <it>Richard P. Brent</it>,
146 * Dover 0-486-41998-3
147 * </blockquote>
148 *
149 * @param lo Lower bound of the search interval.
150 * @param hi Higher bound of the search interval.
151 * @param fLo Function value at the lower bound of the search interval.
152 * @param fHi Function value at the higher bound of the search interval.
153 * @return the value where the function is zero.
154 */
155 private double brent(double lo, double hi,
156 double fLo, double fHi) {
157 double a = lo;
158 double fa = fLo;
159 double b = hi;
160 double fb = fHi;
161 double c = a;
162 double fc = fa;
163 double d = b - a;
164 double e = d;
165
166 final double t = getAbsoluteAccuracy();
167 final double eps = getRelativeAccuracy();
168
169 while (true) {
170 if (FastMath.abs(fc) < FastMath.abs(fb)) {
171 a = b;
172 b = c;
173 c = a;
174 fa = fb;
175 fb = fc;
176 fc = fa;
177 }
178
179 final double tol = 2 * eps * FastMath.abs(b) + t;
180 final double m = 0.5 * (c - b);
181
182 if (FastMath.abs(m) <= tol ||
183 Precision.equals(fb, 0)) {
184 return b;
185 }
186 if (FastMath.abs(e) < tol ||
187 FastMath.abs(fa) <= FastMath.abs(fb)) {
188 // Force bisection.
189 d = m;
190 e = d;
191 } else {
192 double s = fb / fa;
193 double p;
194 double q;
195 // The equality test (a == c) is intentional,
196 // it is part of the original Brent's method and
197 // it should NOT be replaced by proximity test.
198 if (a == c) {
199 // Linear interpolation.
200 p = 2 * m * s;
201 q = 1 - s;
202 } else {
203 // Inverse quadratic interpolation.
204 q = fa / fc;
205 final double r = fb / fc;
206 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
207 q = (q - 1) * (r - 1) * (s - 1);
208 }
209 if (p > 0) {
210 q = -q;
211 } else {
212 p = -p;
213 }
214 s = e;
215 e = d;
216 if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
217 p >= FastMath.abs(0.5 * s * q)) {
218 // Inverse quadratic interpolation gives a value
219 // in the wrong direction, or progress is slow.
220 // Fall back to bisection.
221 d = m;
222 e = d;
223 } else {
224 d = p / q;
225 }
226 }
227 a = b;
228 fa = fb;
229
230 if (FastMath.abs(d) > tol) {
231 b += d;
232 } else if (m > 0) {
233 b += tol;
234 } else {
235 b -= tol;
236 }
237 fb = computeObjectiveValue(b);
238 if ((fb > 0 && fc > 0) ||
239 (fb <= 0 && fc <= 0)) {
240 c = a;
241 fc = fa;
242 d = b - a;
243 e = d;
244 }
245 }
246 }
247 }