Coverage Report - org.apache.commons.math.analysis.MullerSolver

Classes in this File	Line Coverage	Branch Coverage	Complexity
MullerSolver	83%	85%	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		}