Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
MullerSolver |
|
| 9.25;9.25 |
1 | /* |
|
2 | * Copyright 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 | import org.apache.commons.math.ConvergenceException; |
|
19 | import org.apache.commons.math.FunctionEvaluationException; |
|
20 | import org.apache.commons.math.util.MathUtils; |
|
21 | ||
22 | /** |
|
23 | * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html"> |
|
24 | * Muller's Method</a> for root finding of real univariate functions. For |
|
25 | * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, |
|
26 | * chapter 3. |
|
27 | * <p> |
|
28 | * Muller's method applies to both real and complex functions, but here we |
|
29 | * restrict ourselves to real functions. Methods solve() and solve2() find |
|
30 | * real zeros, using different ways to bypass complex arithmetics. |
|
31 | * |
|
32 | * @version $Revision$ $Date$ |
|
33 | */ |
|
34 | public class MullerSolver extends UnivariateRealSolverImpl { |
|
35 | ||
36 | /** serializable version identifier */ |
|
37 | static final long serialVersionUID = 2619993603551148137L; |
|
38 | ||
39 | /** |
|
40 | * Construct a solver for the given function. |
|
41 | * |
|
42 | * @param f function to solve |
|
43 | */ |
|
44 | public MullerSolver(UnivariateRealFunction f) { |
|
45 | 14 | super(f, 100, 1E-6); |
46 | 14 | } |
47 | ||
48 | /** |
|
49 | * Find a real root in the given interval with initial value. |
|
50 | * <p> |
|
51 | * Requires bracketing condition. |
|
52 | * |
|
53 | * @param min the lower bound for the interval |
|
54 | * @param max the upper bound for the interval |
|
55 | * @param initial the start value to use |
|
56 | * @return the point at which the function value is zero |
|
57 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
58 | * or the solver detects convergence problems otherwise |
|
59 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
60 | * function |
|
61 | * @throws IllegalArgumentException if any parameters are invalid |
|
62 | */ |
|
63 | public double solve(double min, double max, double initial) throws |
|
64 | ConvergenceException, FunctionEvaluationException { |
|
65 | ||
66 | // check for zeros before verifying bracketing |
|
67 | 0 | if (f.value(min) == 0.0) { return min; } |
68 | 0 | if (f.value(max) == 0.0) { return max; } |
69 | 0 | if (f.value(initial) == 0.0) { return initial; } |
70 | ||
71 | 0 | verifyBracketing(min, max, f); |
72 | 0 | verifySequence(min, initial, max); |
73 | 0 | if (isBracketing(min, initial, f)) { |
74 | 0 | return solve(min, initial); |
75 | } else { |
|
76 | 0 | return solve(initial, max); |
77 | } |
|
78 | } |
|
79 | ||
80 | /** |
|
81 | * Find a real root in the given interval. |
|
82 | * <p> |
|
83 | * Original Muller's method would have function evaluation at complex point. |
|
84 | * Since our f(x) is real, we have to find ways to avoid that. Bracketing |
|
85 | * condition is one way to go: by requiring bracketing in every iteration, |
|
86 | * the newly computed approximation is guaranteed to be real. |
|
87 | * <p> |
|
88 | * Normally Muller's method converges quadratically in the vicinity of a |
|
89 | * zero, however it may be very slow in regions far away from zeros. For |
|
90 | * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use |
|
91 | * bisection as a safety backup if it performs very poorly. |
|
92 | * <p> |
|
93 | * The formulas here use divided differences directly. |
|
94 | * |
|
95 | * @param min the lower bound for the interval |
|
96 | * @param max the upper bound for the interval |
|
97 | * @return the point at which the function value is zero |
|
98 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
99 | * or the solver detects convergence problems otherwise |
|
100 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
101 | * function |
|
102 | * @throws IllegalArgumentException if any parameters are invalid |
|
103 | */ |
|
104 | public double solve(double min, double max) throws ConvergenceException, |
|
105 | FunctionEvaluationException { |
|
106 | ||
107 | // [x0, x2] is the bracketing interval in each iteration |
|
108 | // x1 is the last approximation and an interpolation point in (x0, x2) |
|
109 | // x is the new root approximation and new x1 for next round |
|
110 | // d01, d12, d012 are divided differences |
|
111 | double x0, x1, x2, x, oldx, y0, y1, y2, y; |
|
112 | double d01, d12, d012, c1, delta, xplus, xminus, tolerance; |
|
113 | ||
114 | 20 | x0 = min; y0 = f.value(x0); |
115 | 20 | x2 = max; y2 = f.value(x2); |
116 | 20 | x1 = 0.5 * (x0 + x2); y1 = f.value(x1); |
117 | ||
118 | // check for zeros before verifying bracketing |
|
119 | 20 | if (y0 == 0.0) { return min; } |
120 | 20 | if (y2 == 0.0) { return max; } |
121 | 20 | verifyBracketing(min, max, f); |
122 | ||
123 | 16 | int i = 1; |
124 | 16 | oldx = Double.POSITIVE_INFINITY; |
125 | 114 | while (i <= maximalIterationCount) { |
126 | // Muller's method employs quadratic interpolation through |
|
127 | // x0, x1, x2 and x is the zero of the interpolating parabola. |
|
128 | // Due to bracketing condition, this parabola must have two |
|
129 | // real roots and we choose one in [x0, x2] to be x. |
|
130 | 114 | d01 = (y1 - y0) / (x1 - x0); |
131 | 114 | d12 = (y2 - y1) / (x2 - x1); |
132 | 114 | d012 = (d12 - d01) / (x2 - x0); |
133 | 114 | c1 = d01 + (x1 - x0) * d012; |
134 | 114 | delta = c1 * c1 - 4 * y1 * d012; |
135 | 114 | xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta)); |
136 | 114 | xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta)); |
137 | // xplus and xminus are two roots of parabola and at least |
|
138 | // one of them should lie in (x0, x2) |
|
139 | 114 | x = isSequence(x0, xplus, x2) ? xplus : xminus; |
140 | 114 | y = f.value(x); |
141 | ||
142 | // check for convergence |
|
143 | 114 | tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); |
144 | 114 | if (Math.abs(x - oldx) <= tolerance) { |
145 | 16 | setResult(x, i); |
146 | 16 | return result; |
147 | } |
|
148 | 98 | if (Math.abs(y) <= functionValueAccuracy) { |
149 | 0 | setResult(x, i); |
150 | 0 | return result; |
151 | } |
|
152 | ||
153 | // Bisect if convergence is too slow. Bisection would waste |
|
154 | // our calculation of x, hopefully it won't happen often. |
|
155 | 98 | boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || |
156 | (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || |
|
157 | (x == x1); |
|
158 | // prepare the new bracketing interval for next iteration |
|
159 | 98 | if (!bisect) { |
160 | 82 | x0 = x < x1 ? x0 : x1; |
161 | 82 | y0 = x < x1 ? y0 : y1; |
162 | 82 | x2 = x > x1 ? x2 : x1; |
163 | 82 | y2 = x > x1 ? y2 : y1; |
164 | 82 | x1 = x; y1 = y; |
165 | 82 | oldx = x; |
166 | } else { |
|
167 | 16 | double xm = 0.5 * (x0 + x2); |
168 | 16 | double ym = f.value(xm); |
169 | 16 | if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) { |
170 | 12 | x2 = xm; y2 = ym; |
171 | } else { |
|
172 | 4 | x0 = xm; y0 = ym; |
173 | } |
|
174 | 16 | x1 = 0.5 * (x0 + x2); |
175 | 16 | y1 = f.value(x1); |
176 | 16 | oldx = Double.POSITIVE_INFINITY; |
177 | } |
|
178 | 98 | i++; |
179 | } |
|
180 | 0 | throw new ConvergenceException("Maximum number of iterations exceeded."); |
181 | } |
|
182 | ||
183 | /** |
|
184 | * Find a real root in the given interval. |
|
185 | * <p> |
|
186 | * solve2() differs from solve() in the way it avoids complex operations. |
|
187 | * Except for the initial [min, max], solve2() does not require bracketing |
|
188 | * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex |
|
189 | * number arises in the computation, we simply use its modulus as real |
|
190 | * approximation. |
|
191 | * <p> |
|
192 | * Because the interval may not be bracketing, bisection alternative is |
|
193 | * not applicable here. However in practice our treatment usually works |
|
194 | * well, especially near real zeros where the imaginary part of complex |
|
195 | * approximation is often negligible. |
|
196 | * <p> |
|
197 | * The formulas here do not use divided differences directly. |
|
198 | * |
|
199 | * @param min the lower bound for the interval |
|
200 | * @param max the upper bound for the interval |
|
201 | * @return the point at which the function value is zero |
|
202 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
203 | * or the solver detects convergence problems otherwise |
|
204 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
205 | * function |
|
206 | * @throws IllegalArgumentException if any parameters are invalid |
|
207 | */ |
|
208 | public double solve2(double min, double max) throws ConvergenceException, |
|
209 | FunctionEvaluationException { |
|
210 | ||
211 | // x2 is the last root approximation |
|
212 | // x is the new approximation and new x2 for next round |
|
213 | // x0 < x1 < x2 does not hold here |
|
214 | double x0, x1, x2, x, oldx, y0, y1, y2, y; |
|
215 | double q, A, B, C, delta, denominator, tolerance; |
|
216 | ||
217 | 16 | x0 = min; y0 = f.value(x0); |
218 | 16 | x1 = max; y1 = f.value(x1); |
219 | 16 | x2 = 0.5 * (x0 + x1); y2 = f.value(x2); |
220 | ||
221 | // check for zeros before verifying bracketing |
|
222 | 16 | if (y0 == 0.0) { return min; } |
223 | 16 | if (y1 == 0.0) { return max; } |
224 | 16 | verifyBracketing(min, max, f); |
225 | ||
226 | 16 | int i = 1; |
227 | 16 | oldx = Double.POSITIVE_INFINITY; |
228 | 218 | while (i <= maximalIterationCount) { |
229 | // quadratic interpolation through x0, x1, x2 |
|
230 | 218 | q = (x2 - x1) / (x1 - x0); |
231 | 218 | A = q * (y2 - (1 + q) * y1 + q * y0); |
232 | 218 | B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; |
233 | 218 | C = (1 + q) * y2; |
234 | 218 | delta = B * B - 4 * A * C; |
235 | 218 | if (delta >= 0.0) { |
236 | // choose a denominator larger in magnitude |
|
237 | 128 | double dplus = B + Math.sqrt(delta); |
238 | 128 | double dminus = B - Math.sqrt(delta); |
239 | 128 | denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus; |
240 | } else { |
|
241 | // take the modulus of (B +/- Math.sqrt(delta)) |
|
242 | 90 | denominator = Math.sqrt(B * B - delta); |
243 | } |
|
244 | 218 | if (denominator != 0) { |
245 | 218 | x = x2 - 2.0 * C * (x2 - x1) / denominator; |
246 | // perturb x if it coincides with x1 or x2 |
|
247 | 222 | while (x == x1 || x == x2) { |
248 | 4 | x += absoluteAccuracy; |
249 | } |
|
250 | } else { |
|
251 | // extremely rare case, get a random number to skip it |
|
252 | 0 | x = min + Math.random() * (max - min); |
253 | 0 | oldx = Double.POSITIVE_INFINITY; |
254 | } |
|
255 | 218 | y = f.value(x); |
256 | ||
257 | // check for convergence |
|
258 | 218 | tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); |
259 | 218 | if (Math.abs(x - oldx) <= tolerance) { |
260 | 16 | setResult(x, i); |
261 | 16 | return result; |
262 | } |
|
263 | 202 | if (Math.abs(y) <= functionValueAccuracy) { |
264 | 0 | setResult(x, i); |
265 | 0 | return result; |
266 | } |
|
267 | ||
268 | // prepare the next iteration |
|
269 | 202 | x0 = x1; y0 = y1; |
270 | 202 | x1 = x2; y1 = y2; |
271 | 202 | x2 = x; y2 = y; |
272 | 202 | oldx = x; |
273 | 202 | i++; |
274 | } |
|
275 | 0 | throw new ConvergenceException("Maximum number of iterations exceeded."); |
276 | } |
|
277 | } |