1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.analysis.solvers;
18
19
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.MathRuntimeException;
22 import org.apache.commons.math.MaxIterationsExceededException;
23 import org.apache.commons.math.analysis.UnivariateRealFunction;
24
25
26
27
28
29
30
31
32
33 public class BrentSolver extends UnivariateRealSolverImpl {
34
35
36 private static final long serialVersionUID = 7694577816772532779L;
37
38
39
40
41
42
43
44
45
46
47 @Deprecated
48 public BrentSolver(UnivariateRealFunction f) {
49 super(f, 100, 1E-6);
50 }
51
52
53
54
55 public BrentSolver() {
56 super(100, 1E-6);
57 }
58
59
60 @Deprecated
61 public double solve(double min, double max)
62 throws MaxIterationsExceededException, FunctionEvaluationException {
63 return solve(f, min, max);
64 }
65
66
67 @Deprecated
68 public double solve(double min, double max, double initial)
69 throws MaxIterationsExceededException, FunctionEvaluationException {
70 return solve(f, min, max, initial);
71 }
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93 public double solve(final UnivariateRealFunction f,
94 final double min, final double max, final double initial)
95 throws MaxIterationsExceededException, FunctionEvaluationException {
96
97 clearResult();
98 verifySequence(min, initial, max);
99
100
101 double yInitial = f.value(initial);
102 if (Math.abs(yInitial) <= functionValueAccuracy) {
103 setResult(initial, 0);
104 return result;
105 }
106
107
108 double yMin = f.value(min);
109 if (Math.abs(yMin) <= functionValueAccuracy) {
110 setResult(yMin, 0);
111 return result;
112 }
113
114
115 if (yInitial * yMin < 0) {
116 return solve(f, min, yMin, initial, yInitial, min, yMin);
117 }
118
119
120 double yMax = f.value(max);
121 if (Math.abs(yMax) <= functionValueAccuracy) {
122 setResult(yMax, 0);
123 return result;
124 }
125
126
127 if (yInitial * yMax < 0) {
128 return solve(f, initial, yInitial, max, yMax, initial, yInitial);
129 }
130
131
132 return solve(f, min, yMin, max, yMax, initial, yInitial);
133
134 }
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153 public double solve(final UnivariateRealFunction f,
154 final double min, final double max)
155 throws MaxIterationsExceededException,
156 FunctionEvaluationException {
157
158 clearResult();
159 verifyInterval(min, max);
160
161 double ret = Double.NaN;
162
163 double yMin = f.value(min);
164 double yMax = f.value(max);
165
166
167 double sign = yMin * yMax;
168 if (sign > 0) {
169
170 if (Math.abs(yMin) <= functionValueAccuracy) {
171 setResult(min, 0);
172 ret = min;
173 } else if (Math.abs(yMax) <= functionValueAccuracy) {
174 setResult(max, 0);
175 ret = max;
176 } else {
177
178 throw MathRuntimeException.createIllegalArgumentException(
179 "function values at endpoints do not have different signs. " +
180 "Endpoints: [{0}, {1}], Values: [{2}, {3}]",
181 min, max, yMin, yMax);
182 }
183 } else if (sign < 0){
184
185 ret = solve(f, min, yMin, max, yMax, min, yMin);
186 } else {
187
188 if (yMin == 0.0) {
189 ret = min;
190 } else {
191 ret = max;
192 }
193 }
194
195 return ret;
196 }
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215 private double solve(final UnivariateRealFunction f,
216 double x0, double y0,
217 double x1, double y1,
218 double x2, double y2)
219 throws MaxIterationsExceededException, FunctionEvaluationException {
220
221 double delta = x1 - x0;
222 double oldDelta = delta;
223
224 int i = 0;
225 while (i < maximalIterationCount) {
226 if (Math.abs(y2) < Math.abs(y1)) {
227
228 x0 = x1;
229 x1 = x2;
230 x2 = x0;
231 y0 = y1;
232 y1 = y2;
233 y2 = y0;
234 }
235 if (Math.abs(y1) <= functionValueAccuracy) {
236
237
238
239 setResult(x1, i);
240 return result;
241 }
242 double dx = (x2 - x1);
243 double tolerance =
244 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
245 if (Math.abs(dx) <= tolerance) {
246 setResult(x1, i);
247 return result;
248 }
249 if ((Math.abs(oldDelta) < tolerance) ||
250 (Math.abs(y0) <= Math.abs(y1))) {
251
252 delta = 0.5 * dx;
253 oldDelta = delta;
254 } else {
255 double r3 = y1 / y0;
256 double p;
257 double p1;
258
259
260
261 if (x0 == x2) {
262
263 p = dx * r3;
264 p1 = 1.0 - r3;
265 } else {
266
267 double r1 = y0 / y2;
268 double r2 = y1 / y2;
269 p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
270 p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
271 }
272 if (p > 0.0) {
273 p1 = -p1;
274 } else {
275 p = -p;
276 }
277 if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
278 p >= Math.abs(0.5 * oldDelta * p1)) {
279
280
281
282 delta = 0.5 * dx;
283 oldDelta = delta;
284 } else {
285 oldDelta = delta;
286 delta = p / p1;
287 }
288 }
289
290 x0 = x1;
291 y0 = y1;
292
293 if (Math.abs(delta) > tolerance) {
294 x1 = x1 + delta;
295 } else if (dx > 0.0) {
296 x1 = x1 + 0.5 * tolerance;
297 } else if (dx <= 0.0) {
298 x1 = x1 - 0.5 * tolerance;
299 }
300 y1 = f.value(x1);
301 if ((y1 > 0) == (y2 > 0)) {
302 x2 = x0;
303 y2 = y0;
304 delta = x1 - x0;
305 oldDelta = delta;
306 }
307 i++;
308 }
309 throw new MaxIterationsExceededException(maximalIterationCount);
310 }
311 }