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
9
* http://www.apache.org/licenses/LICENSE-2.0
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.
17
package org.apache.commons.math.analysis.solvers;
19
import org.apache.commons.math.ConvergenceException;
20
import org.apache.commons.math.FunctionEvaluationException;
21
import org.apache.commons.math.MaxIterationsExceededException;
22
import org.apache.commons.math.analysis.UnivariateRealFunction;
23
import org.apache.commons.math.util.MathUtils;
26
* Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
27
* Muller's Method</a> for root finding of real univariate functions. For
28
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
31
* Muller's method applies to both real and complex functions, but here we
32
* restrict ourselves to real functions. Methods solve() and solve2() find
33
* real zeros, using different ways to bypass complex arithmetics.</p>
35
* @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
38
public class MullerSolver extends UnivariateRealSolverImpl {
41
* Construct a solver for the given function.
43
* @param f function to solve
44
* @deprecated as of 2.0 the function to solve is passed as an argument
45
* to the {@link #solve(UnivariateRealFunction, double, double)} or
46
* {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
50
public MullerSolver(UnivariateRealFunction f) {
57
public MullerSolver() {
63
public double solve(final double min, final double max)
64
throws ConvergenceException, FunctionEvaluationException {
65
return solve(f, min, max);
70
public double solve(final double min, final double max, final double initial)
71
throws ConvergenceException, FunctionEvaluationException {
72
return solve(f, min, max, initial);
76
* Find a real root in the given interval with initial value.
78
* Requires bracketing condition.</p>
80
* @param f the function to solve
81
* @param min the lower bound for the interval
82
* @param max the upper bound for the interval
83
* @param initial the start value to use
84
* @return the point at which the function value is zero
85
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
86
* or the solver detects convergence problems otherwise
87
* @throws FunctionEvaluationException if an error occurs evaluating the
89
* @throws IllegalArgumentException if any parameters are invalid
91
public double solve(final UnivariateRealFunction f,
92
final double min, final double max, final double initial)
93
throws MaxIterationsExceededException, FunctionEvaluationException {
95
// check for zeros before verifying bracketing
96
if (f.value(min) == 0.0) { return min; }
97
if (f.value(max) == 0.0) { return max; }
98
if (f.value(initial) == 0.0) { return initial; }
100
verifyBracketing(min, max, f);
101
verifySequence(min, initial, max);
102
if (isBracketing(min, initial, f)) {
103
return solve(f, min, initial);
105
return solve(f, initial, max);
110
* Find a real root in the given interval.
112
* Original Muller's method would have function evaluation at complex point.
113
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
114
* condition is one way to go: by requiring bracketing in every iteration,
115
* the newly computed approximation is guaranteed to be real.</p>
117
* Normally Muller's method converges quadratically in the vicinity of a
118
* zero, however it may be very slow in regions far away from zeros. For
119
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
120
* bisection as a safety backup if it performs very poorly.</p>
122
* The formulas here use divided differences directly.</p>
124
* @param f the function to solve
125
* @param min the lower bound for the interval
126
* @param max the upper bound for the interval
127
* @return the point at which the function value is zero
128
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
129
* or the solver detects convergence problems otherwise
130
* @throws FunctionEvaluationException if an error occurs evaluating the
132
* @throws IllegalArgumentException if any parameters are invalid
134
public double solve(final UnivariateRealFunction f,
135
final double min, final double max)
136
throws MaxIterationsExceededException, FunctionEvaluationException {
138
// [x0, x2] is the bracketing interval in each iteration
139
// x1 is the last approximation and an interpolation point in (x0, x2)
140
// x is the new root approximation and new x1 for next round
141
// d01, d12, d012 are divided differences
142
double x0, x1, x2, x, oldx, y0, y1, y2, y;
143
double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
145
x0 = min; y0 = f.value(x0);
146
x2 = max; y2 = f.value(x2);
147
x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
149
// check for zeros before verifying bracketing
150
if (y0 == 0.0) { return min; }
151
if (y2 == 0.0) { return max; }
152
verifyBracketing(min, max, f);
155
oldx = Double.POSITIVE_INFINITY;
156
while (i <= maximalIterationCount) {
157
// Muller's method employs quadratic interpolation through
158
// x0, x1, x2 and x is the zero of the interpolating parabola.
159
// Due to bracketing condition, this parabola must have two
160
// real roots and we choose one in [x0, x2] to be x.
161
d01 = (y1 - y0) / (x1 - x0);
162
d12 = (y2 - y1) / (x2 - x1);
163
d012 = (d12 - d01) / (x2 - x0);
164
c1 = d01 + (x1 - x0) * d012;
165
delta = c1 * c1 - 4 * y1 * d012;
166
xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
167
xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
168
// xplus and xminus are two roots of parabola and at least
169
// one of them should lie in (x0, x2)
170
x = isSequence(x0, xplus, x2) ? xplus : xminus;
173
// check for convergence
174
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
175
if (Math.abs(x - oldx) <= tolerance) {
179
if (Math.abs(y) <= functionValueAccuracy) {
184
// Bisect if convergence is too slow. Bisection would waste
185
// our calculation of x, hopefully it won't happen often.
186
// the real number equality test x == x1 is intentional and
187
// completes the proximity tests above it
188
boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
189
(x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
191
// prepare the new bracketing interval for next iteration
193
x0 = x < x1 ? x0 : x1;
194
y0 = x < x1 ? y0 : y1;
195
x2 = x > x1 ? x2 : x1;
196
y2 = x > x1 ? y2 : y1;
200
double xm = 0.5 * (x0 + x2);
201
double ym = f.value(xm);
202
if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
207
x1 = 0.5 * (x0 + x2);
209
oldx = Double.POSITIVE_INFINITY;
213
throw new MaxIterationsExceededException(maximalIterationCount);
217
* Find a real root in the given interval.
219
* solve2() differs from solve() in the way it avoids complex operations.
220
* Except for the initial [min, max], solve2() does not require bracketing
221
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
222
* number arises in the computation, we simply use its modulus as real
225
* Because the interval may not be bracketing, bisection alternative is
226
* not applicable here. However in practice our treatment usually works
227
* well, especially near real zeros where the imaginary part of complex
228
* approximation is often negligible.</p>
230
* The formulas here do not use divided differences directly.</p>
232
* @param min the lower bound for the interval
233
* @param max the upper bound for the interval
234
* @return the point at which the function value is zero
235
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
236
* or the solver detects convergence problems otherwise
237
* @throws FunctionEvaluationException if an error occurs evaluating the
239
* @throws IllegalArgumentException if any parameters are invalid
240
* @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)}
244
public double solve2(final double min, final double max)
245
throws MaxIterationsExceededException, FunctionEvaluationException {
246
return solve2(f, min, max);
250
* Find a real root in the given interval.
252
* solve2() differs from solve() in the way it avoids complex operations.
253
* Except for the initial [min, max], solve2() does not require bracketing
254
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
255
* number arises in the computation, we simply use its modulus as real
258
* Because the interval may not be bracketing, bisection alternative is
259
* not applicable here. However in practice our treatment usually works
260
* well, especially near real zeros where the imaginary part of complex
261
* approximation is often negligible.</p>
263
* The formulas here do not use divided differences directly.</p>
265
* @param f the function to solve
266
* @param min the lower bound for the interval
267
* @param max the upper bound for the interval
268
* @return the point at which the function value is zero
269
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
270
* or the solver detects convergence problems otherwise
271
* @throws FunctionEvaluationException if an error occurs evaluating the
273
* @throws IllegalArgumentException if any parameters are invalid
275
public double solve2(final UnivariateRealFunction f,
276
final double min, final double max)
277
throws MaxIterationsExceededException, FunctionEvaluationException {
279
// x2 is the last root approximation
280
// x is the new approximation and new x2 for next round
281
// x0 < x1 < x2 does not hold here
282
double x0, x1, x2, x, oldx, y0, y1, y2, y;
283
double q, A, B, C, delta, denominator, tolerance;
285
x0 = min; y0 = f.value(x0);
286
x1 = max; y1 = f.value(x1);
287
x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
289
// check for zeros before verifying bracketing
290
if (y0 == 0.0) { return min; }
291
if (y1 == 0.0) { return max; }
292
verifyBracketing(min, max, f);
295
oldx = Double.POSITIVE_INFINITY;
296
while (i <= maximalIterationCount) {
297
// quadratic interpolation through x0, x1, x2
298
q = (x2 - x1) / (x1 - x0);
299
A = q * (y2 - (1 + q) * y1 + q * y0);
300
B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
302
delta = B * B - 4 * A * C;
304
// choose a denominator larger in magnitude
305
double dplus = B + Math.sqrt(delta);
306
double dminus = B - Math.sqrt(delta);
307
denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
309
// take the modulus of (B +/- Math.sqrt(delta))
310
denominator = Math.sqrt(B * B - delta);
312
if (denominator != 0) {
313
x = x2 - 2.0 * C * (x2 - x1) / denominator;
314
// perturb x if it exactly coincides with x1 or x2
315
// the equality tests here are intentional
316
while (x == x1 || x == x2) {
317
x += absoluteAccuracy;
320
// extremely rare case, get a random number to skip it
321
x = min + Math.random() * (max - min);
322
oldx = Double.POSITIVE_INFINITY;
326
// check for convergence
327
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
328
if (Math.abs(x - oldx) <= tolerance) {
332
if (Math.abs(y) <= functionValueAccuracy) {
337
// prepare the next iteration
344
throw new MaxIterationsExceededException(maximalIterationCount);