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