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 * http://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 package org.apache.commons.math.analysis.solvers;
18
19 import org.apache.commons.math.util.FastMath;
20
21 /**
22 * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
23 * Muller's Method</a> for root finding of real univariate functions. For
24 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
25 * chapter 3.
26 * <p>
27 * Muller's method applies to both real and complex functions, but here we
28 * restrict ourselves to real functions.
29 * This class differs from {@link MullerSolver} in the way it avoids complex
30 * operations.</p>
31 * Muller's original method would have function evaluation at complex point.
32 * Since our f(x) is real, we have to find ways to avoid that. Bracketing
33 * condition is one way to go: by requiring bracketing in every iteration,
34 * the newly computed approximation is guaranteed to be real.</p>
35 * <p>
36 * Normally Muller's method converges quadratically in the vicinity of a
37 * zero, however it may be very slow in regions far away from zeros. For
38 * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
39 * bisection as a safety backup if it performs very poorly.</p>
40 * <p>
41 * The formulas here use divided differences directly.</p>
42 *
43 * @version $Id: MullerSolver.java 1183138 2011-10-13 22:21:04Z erans $
44 * @since 1.2
45 * @see MullerSolver2
46 */
47 public class MullerSolver extends AbstractUnivariateRealSolver {
48
49 /** Default absolute accuracy. */
50 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
51
52 /**
53 * Construct a solver with default accuracy (1e-6).
54 */
55 public MullerSolver() {
56 this(DEFAULT_ABSOLUTE_ACCURACY);
57 }
58 /**
59 * Construct a solver.
60 *
61 * @param absoluteAccuracy Absolute accuracy.
62 */
63 public MullerSolver(double absoluteAccuracy) {
64 super(absoluteAccuracy);
65 }
66 /**
67 * Construct a solver.
68 *
69 * @param relativeAccuracy Relative accuracy.
70 * @param absoluteAccuracy Absolute accuracy.
71 */
72 public MullerSolver(double relativeAccuracy,
73 double absoluteAccuracy) {
74 super(relativeAccuracy, absoluteAccuracy);
75 }
76
77 /**
78 * {@inheritDoc}
79 */
80 @Override
81 protected double doSolve() {
82 final double min = getMin();
83 final double max = getMax();
84 final double initial = getStartValue();
85
86 final double functionValueAccuracy = getFunctionValueAccuracy();
87
88 verifySequence(min, initial, max);
89
90 // check for zeros before verifying bracketing
91 final double fMin = computeObjectiveValue(min);
92 if (FastMath.abs(fMin) < functionValueAccuracy) {
93 return min;
94 }
95 final double fMax = computeObjectiveValue(max);
96 if (FastMath.abs(fMax) < functionValueAccuracy) {
97 return max;
98 }
99 final double fInitial = computeObjectiveValue(initial);
100 if (FastMath.abs(fInitial) < functionValueAccuracy) {
101 return initial;
102 }
103
104 verifyBracketing(min, max);
105
106 if (isBracketing(min, initial)) {
107 return solve(min, initial, fMin, fInitial);
108 } else {
109 return solve(initial, max, fInitial, fMax);
110 }
111 }
112
113 /**
114 * Find a real root in the given interval.
115 *
116 * @param min Lower bound for the interval.
117 * @param max Upper bound for the interval.
118 * @param fMin function value at the lower bound.
119 * @param fMax function value at the upper bound.
120 * @return the point at which the function value is zero.
121 */
122 private double solve(double min, double max,
123 double fMin, double fMax) {
124 final double relativeAccuracy = getRelativeAccuracy();
125 final double absoluteAccuracy = getAbsoluteAccuracy();
126 final double functionValueAccuracy = getFunctionValueAccuracy();
127
128 // [x0, x2] is the bracketing interval in each iteration
129 // x1 is the last approximation and an interpolation point in (x0, x2)
130 // x is the new root approximation and new x1 for next round
131 // d01, d12, d012 are divided differences
132
133 double x0 = min;
134 double y0 = fMin;
135 double x2 = max;
136 double y2 = fMax;
137 double x1 = 0.5 * (x0 + x2);
138 double y1 = computeObjectiveValue(x1);
139
140 double oldx = Double.POSITIVE_INFINITY;
141 while (true) {
142 // Muller's method employs quadratic interpolation through
143 // x0, x1, x2 and x is the zero of the interpolating parabola.
144 // Due to bracketing condition, this parabola must have two
145 // real roots and we choose one in [x0, x2] to be x.
146 final double d01 = (y1 - y0) / (x1 - x0);
147 final double d12 = (y2 - y1) / (x2 - x1);
148 final double d012 = (d12 - d01) / (x2 - x0);
149 final double c1 = d01 + (x1 - x0) * d012;
150 final double delta = c1 * c1 - 4 * y1 * d012;
151 final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta));
152 final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta));
153 // xplus and xminus are two roots of parabola and at least
154 // one of them should lie in (x0, x2)
155 final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
156 final double y = computeObjectiveValue(x);
157
158 // check for convergence
159 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
160 if (FastMath.abs(x - oldx) <= tolerance ||
161 FastMath.abs(y) <= functionValueAccuracy) {
162 return x;
163 }
164
165 // Bisect if convergence is too slow. Bisection would waste
166 // our calculation of x, hopefully it won't happen often.
167 // the real number equality test x == x1 is intentional and
168 // completes the proximity tests above it
169 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
170 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
171 (x == x1);
172 // prepare the new bracketing interval for next iteration
173 if (!bisect) {
174 x0 = x < x1 ? x0 : x1;
175 y0 = x < x1 ? y0 : y1;
176 x2 = x > x1 ? x2 : x1;
177 y2 = x > x1 ? y2 : y1;
178 x1 = x; y1 = y;
179 oldx = x;
180 } else {
181 double xm = 0.5 * (x0 + x2);
182 double ym = computeObjectiveValue(xm);
183 if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) {
184 x2 = xm; y2 = ym;
185 } else {
186 x0 = xm; y0 = ym;
187 }
188 x1 = 0.5 * (x0 + x2);
189 y1 = computeObjectiveValue(x1);
190 oldx = Double.POSITIVE_INFINITY;
191 }
192 }
193 }
194 }