Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
BrentSolver |
|
| 6.666666666666667;6.667 |
1 | /* |
|
2 | * Copyright 2003-2005 The Apache Software Foundation. |
|
3 | * |
|
4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
|
5 | * you may not use this file except in compliance with the License. |
|
6 | * You may obtain a copy of the License at |
|
7 | * |
|
8 | * http://www.apache.org/licenses/LICENSE-2.0 |
|
9 | * |
|
10 | * Unless required by applicable law or agreed to in writing, software |
|
11 | * distributed under the License is distributed on an "AS IS" BASIS, |
|
12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
13 | * See the License for the specific language governing permissions and |
|
14 | * limitations under the License. |
|
15 | */ |
|
16 | package org.apache.commons.math.analysis; |
|
17 | ||
18 | ||
19 | import org.apache.commons.math.ConvergenceException; |
|
20 | import org.apache.commons.math.FunctionEvaluationException; |
|
21 | ||
22 | /** |
|
23 | * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html"> |
|
24 | * Brent algorithm</a> for finding zeros of real univariate functions. |
|
25 | * <p> |
|
26 | * The function should be continuous but not necessarily smooth. |
|
27 | * |
|
28 | * @version $Revision$ $Date: 2005-06-03 22:36:42 -0700 (Fri, 03 Jun 2005) $ |
|
29 | */ |
|
30 | public class BrentSolver extends UnivariateRealSolverImpl { |
|
31 | ||
32 | /** Serializable version identifier */ |
|
33 | static final long serialVersionUID = 3350616277306882875L; |
|
34 | ||
35 | /** |
|
36 | * Construct a solver for the given function. |
|
37 | * |
|
38 | * @param f function to solve. |
|
39 | */ |
|
40 | public BrentSolver(UnivariateRealFunction f) { |
|
41 | 204 | super(f, 100, 1E-6); |
42 | 202 | } |
43 | ||
44 | /** |
|
45 | * Find a zero in the given interval. |
|
46 | * <p> |
|
47 | * Throws <code>ConvergenceException</code> if the values of the function |
|
48 | * at the endpoints of the interval have the same sign. |
|
49 | * |
|
50 | * @param min the lower bound for the interval. |
|
51 | * @param max the upper bound for the interval. |
|
52 | * @param initial the start value to use (ignored). |
|
53 | * @return the value where the function is zero |
|
54 | * @throws ConvergenceException the maximum iteration count is exceeded |
|
55 | * @throws FunctionEvaluationException if an error occurs evaluating |
|
56 | * the function |
|
57 | * @throws IllegalArgumentException if initial is not between min and max |
|
58 | */ |
|
59 | public double solve(double min, double max, double initial) |
|
60 | throws ConvergenceException, FunctionEvaluationException { |
|
61 | ||
62 | 0 | return solve(min, max); |
63 | } |
|
64 | ||
65 | /** |
|
66 | * Find a zero in the given interval. |
|
67 | * <p> |
|
68 | * Requires that the values of the function at the endpoints have opposite |
|
69 | * signs. An <code>IllegalArgumentException</code> is thrown if this is not |
|
70 | * the case. |
|
71 | * |
|
72 | * @param min the lower bound for the interval. |
|
73 | * @param max the upper bound for the interval. |
|
74 | * @return the value where the function is zero |
|
75 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
76 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
77 | * function |
|
78 | * @throws IllegalArgumentException if min is not less than max or the |
|
79 | * signs of the values of the function at the endpoints are not opposites |
|
80 | */ |
|
81 | public double solve(double min, double max) throws ConvergenceException, |
|
82 | FunctionEvaluationException { |
|
83 | ||
84 | 224 | clearResult(); |
85 | 224 | verifyInterval(min, max); |
86 | ||
87 | // Index 0 is the old approximation for the root. |
|
88 | // Index 1 is the last calculated approximation for the root. |
|
89 | // Index 2 is a bracket for the root with respect to x1. |
|
90 | 222 | double x0 = min; |
91 | 222 | double x1 = max; |
92 | double y0; |
|
93 | double y1; |
|
94 | 222 | y0 = f.value(x0); |
95 | 222 | y1 = f.value(x1); |
96 | ||
97 | // Verify bracketing |
|
98 | 222 | if (y0 * y1 >= 0) { |
99 | 8 | throw new IllegalArgumentException |
100 | ("Function values at endpoints do not have different signs." + |
|
101 | " Endpoints: [" + min + "," + max + "]" + |
|
102 | " Values: [" + y0 + "," + y1 + "]"); |
|
103 | } |
|
104 | ||
105 | 214 | double x2 = x0; |
106 | 214 | double y2 = y0; |
107 | 214 | double delta = x1 - x0; |
108 | 214 | double oldDelta = delta; |
109 | ||
110 | 214 | int i = 0; |
111 | 1602 | while (i < maximalIterationCount) { |
112 | 1602 | if (Math.abs(y2) < Math.abs(y1)) { |
113 | 404 | x0 = x1; |
114 | 404 | x1 = x2; |
115 | 404 | x2 = x0; |
116 | 404 | y0 = y1; |
117 | 404 | y1 = y2; |
118 | 404 | y2 = y0; |
119 | } |
|
120 | 1602 | if (Math.abs(y1) <= functionValueAccuracy) { |
121 | // Avoid division by very small values. Assume |
|
122 | // the iteration has converged (the problem may |
|
123 | // still be ill conditioned) |
|
124 | 12 | setResult(x1, i); |
125 | 12 | return result; |
126 | } |
|
127 | 1590 | double dx = (x2 - x1); |
128 | 1590 | double tolerance = |
129 | Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy); |
|
130 | 1590 | if (Math.abs(dx) <= tolerance) { |
131 | 202 | setResult(x1, i); |
132 | 202 | return result; |
133 | } |
|
134 | 1388 | if ((Math.abs(oldDelta) < tolerance) || |
135 | (Math.abs(y0) <= Math.abs(y1))) { |
|
136 | // Force bisection. |
|
137 | 18 | delta = 0.5 * dx; |
138 | 18 | oldDelta = delta; |
139 | } else { |
|
140 | 1370 | double r3 = y1 / y0; |
141 | double p; |
|
142 | double p1; |
|
143 | 1370 | if (x0 == x2) { |
144 | // Linear interpolation. |
|
145 | 924 | p = dx * r3; |
146 | 924 | p1 = 1.0 - r3; |
147 | } else { |
|
148 | // Inverse quadratic interpolation. |
|
149 | 446 | double r1 = y0 / y2; |
150 | 446 | double r2 = y1 / y2; |
151 | 446 | p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0)); |
152 | 446 | p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0); |
153 | } |
|
154 | 1370 | if (p > 0.0) { |
155 | 670 | p1 = -p1; |
156 | } else { |
|
157 | 700 | p = -p; |
158 | } |
|
159 | 1370 | if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) || |
160 | p >= Math.abs(0.5 * oldDelta * p1)) { |
|
161 | // Inverse quadratic interpolation gives a value |
|
162 | // in the wrong direction, or progress is slow. |
|
163 | // Fall back to bisection. |
|
164 | 64 | delta = 0.5 * dx; |
165 | 64 | oldDelta = delta; |
166 | } else { |
|
167 | 1306 | oldDelta = delta; |
168 | 1306 | delta = p / p1; |
169 | } |
|
170 | } |
|
171 | // Save old X1, Y1 |
|
172 | 1388 | x0 = x1; |
173 | 1388 | y0 = y1; |
174 | // Compute new X1, Y1 |
|
175 | 1388 | if (Math.abs(delta) > tolerance) { |
176 | 1172 | x1 = x1 + delta; |
177 | 216 | } else if (dx > 0.0) { |
178 | 122 | x1 = x1 + 0.5 * tolerance; |
179 | 94 | } else if (dx <= 0.0) { |
180 | 94 | x1 = x1 - 0.5 * tolerance; |
181 | } |
|
182 | 1388 | y1 = f.value(x1); |
183 | 1388 | if ((y1 > 0) == (y2 > 0)) { |
184 | 920 | x2 = x0; |
185 | 920 | y2 = y0; |
186 | 920 | delta = x1 - x0; |
187 | 920 | oldDelta = delta; |
188 | } |
|
189 | 1388 | i++; |
190 | } |
|
191 | 0 | throw new ConvergenceException("Maximum number of iterations exceeded."); |
192 | } |
|
193 | } |