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

Classes in this File	Line Coverage	Branch Coverage	Complexity
LaguerreSolver	87%	83%	5.625;5.625

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.complex.*;
21
22		/**
23		* Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
24		* Laguerre's Method</a> for root finding of real coefficient polynomials.
25		* For reference, see <b>A First Course in Numerical Analysis</b>,
26		* ISBN 048641454X, chapter 8.
27		* <p>
28		* Laguerre's method is global in the sense that it can start with any initial
29		* approximation and be able to solve all roots from that point.
30		*
31		* @version $Revision$ $Date$
32		*/
33		public class LaguerreSolver extends UnivariateRealSolverImpl {
34
35		/** serializable version identifier */
36		static final long serialVersionUID = 5287689975005870178L;
37
38		/** polynomial function to solve */
39		private PolynomialFunction p;
40
41		/**
42		* Construct a solver for the given function.
43		*
44		* @param f function to solve
45		* @throws IllegalArgumentException if function is not polynomial
46		*/
47		public LaguerreSolver(UnivariateRealFunction f) throws
48		IllegalArgumentException {
49
50	12	super(f, 100, 1E-6);
51	12	if (f instanceof PolynomialFunction) {
52	10	p = (PolynomialFunction)f;
53		} else {
54	2	throw new IllegalArgumentException("Function is not polynomial.");
55		}
56	10	}
57
58		/**
59		* Returns a copy of the polynomial function.
60		*
61		* @return a fresh copy of the polynomial function
62		*/
63		public PolynomialFunction getPolynomialFunction() {
64	0	return new PolynomialFunction(p.getCoefficients());
65		}
66
67		/**
68		* Find a real root in the given interval with initial value.
69		* <p>
70		* Requires bracketing condition.
71		*
72		* @param min the lower bound for the interval
73		* @param max the upper bound for the interval
74		* @param initial the start value to use
75		* @return the point at which the function value is zero
76		* @throws ConvergenceException if the maximum iteration count is exceeded
77		* or the solver detects convergence problems otherwise
78		* @throws FunctionEvaluationException if an error occurs evaluating the
79		* function
80		* @throws IllegalArgumentException if any parameters are invalid
81		*/
82		public double solve(double min, double max, double initial) throws
83		ConvergenceException, FunctionEvaluationException {
84
85		// check for zeros before verifying bracketing
86	0	if (p.value(min) == 0.0) { return min; }
87	0	if (p.value(max) == 0.0) { return max; }
88	0	if (p.value(initial) == 0.0) { return initial; }
89
90	0	verifyBracketing(min, max, p);
91	0	verifySequence(min, initial, max);
92	0	if (isBracketing(min, initial, p)) {
93	0	return solve(min, initial);
94		} else {
95	0	return solve(initial, max);
96		}
97		}
98
99		/**
100		* Find a real root in the given interval.
101		* <p>
102		* Despite the bracketing condition, the root returned by solve(Complex[],
103		* Complex) may not be a real zero inside [min, max]. For example,
104		* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
105		* another initial value, or, as we did here, call solveAll() to obtain
106		* all roots and pick up the one that we're looking for.
107		*
108		* @param min the lower bound for the interval
109		* @param max the upper bound for the interval
110		* @return the point at which the function value is zero
111		* @throws ConvergenceException if the maximum iteration count is exceeded
112		* or the solver detects convergence problems otherwise
113		* @throws FunctionEvaluationException if an error occurs evaluating the
114		* function
115		* @throws IllegalArgumentException if any parameters are invalid
116		*/
117		public double solve(double min, double max) throws ConvergenceException,
118		FunctionEvaluationException {
119
120		Complex c[], root[], initial, z;
121
122		// check for zeros before verifying bracketing
123	16	if (p.value(min) == 0.0) { return min; }
124	16	if (p.value(max) == 0.0) { return max; }
125	16	verifyBracketing(min, max, p);
126
127	12	double coefficients[] = p.getCoefficients();
128	12	c = new Complex[coefficients.length];
129	64	for (int i = 0; i < coefficients.length; i++) {
130	52	c[i] = new Complex(coefficients[i], 0.0);
131		}
132	12	initial = new Complex(0.5 * (min + max), 0.0);
133	12	z = solve(c, initial);
134	12	if (isRootOK(min, max, z)) {
135	10	setResult(z.getReal(), iterationCount);
136	10	return result;
137		}
138
139		// solve all roots and select the one we're seeking
140	2	root = solveAll(c, initial);
141	10	for (int i = 0; i < root.length; i++) {
142	10	if (isRootOK(min, max, root[i])) {
143	2	setResult(root[i].getReal(), iterationCount);
144	2	return result;
145		}
146		}
147
148		// should never happen
149	0	throw new ConvergenceException("Convergence failed.");
150		}
151
152		/**
153		* Returns true iff the given complex root is actually a real zero
154		* in the given interval, within the solver tolerance level.
155		*
156		* @param min the lower bound for the interval
157		* @param max the upper bound for the interval
158		* @param z the complex root
159		* @return true iff z is the sought-after real zero
160		*/
161		protected boolean isRootOK(double min, double max, Complex z) {
162	22	double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
163	22	return (isSequence(min, z.getReal(), max)) &&
164		(Math.abs(z.getImaginary()) <= tolerance ||
165		z.abs() <= functionValueAccuracy);
166		}
167
168		/**
169		* Find all complex roots for the polynomial with the given coefficients,
170		* starting from the given initial value.
171		*
172		* @param coefficients the polynomial coefficients array
173		* @param initial the start value to use
174		* @return the point at which the function value is zero
175		* @throws ConvergenceException if the maximum iteration count is exceeded
176		* or the solver detects convergence problems otherwise
177		* @throws FunctionEvaluationException if an error occurs evaluating the
178		* function
179		* @throws IllegalArgumentException if any parameters are invalid
180		*/
181		public Complex[] solveAll(double coefficients[], double initial) throws
182		ConvergenceException, FunctionEvaluationException {
183
184	2	Complex c[] = new Complex[coefficients.length];
185	2	Complex z = new Complex(initial, 0.0);
186	14	for (int i = 0; i < c.length; i++) {
187	12	c[i] = new Complex(coefficients[i], 0.0);
188		}
189	2	return solveAll(c, z);
190		}
191
192		/**
193		* Find all complex roots for the polynomial with the given coefficients,
194		* starting from the given initial value.
195		*
196		* @param coefficients the polynomial coefficients array
197		* @param initial the start value to use
198		* @return the point at which the function value is zero
199		* @throws ConvergenceException if the maximum iteration count is exceeded
200		* or the solver detects convergence problems otherwise
201		* @throws FunctionEvaluationException if an error occurs evaluating the
202		* function
203		* @throws IllegalArgumentException if any parameters are invalid
204		*/
205		public Complex[] solveAll(Complex coefficients[], Complex initial) throws
206		ConvergenceException, FunctionEvaluationException {
207
208	4	int i, j, n, iterationCount = 0;
209		Complex root[], c[], subarray[], oldc, newc;
210
211	4	n = coefficients.length - 1;
212	4	if (n < 1) {
213	0	throw new IllegalArgumentException
214		("Polynomial degree must be positive: degree=" + n);
215		}
216	4	c = new Complex[n+1]; // coefficients for deflated polynomial
217	28	for (i = 0; i <= n; i++) {
218	24	c[i] = coefficients[i];
219		}
220
221		// solve individual root successively
222	4	root = new Complex[n];
223	24	for (i = 0; i < n; i++) {
224	20	subarray = new Complex[n-i+1];
225	20	System.arraycopy(c, 0, subarray, 0, subarray.length);
226	20	root[i] = solve(subarray, initial);
227		// polynomial deflation using synthetic division
228	20	newc = c[n-i];
229	80	for (j = n-i-1; j >= 0; j--) {
230	60	oldc = c[j];
231	60	c[j] = newc;
232	60	newc = oldc.add(newc.multiply(root[i]));
233		}
234	20	iterationCount += this.iterationCount;
235		}
236
237	4	resultComputed = true;
238	4	this.iterationCount = iterationCount;
239	4	return root;
240		}
241
242		/**
243		* Find a complex root for the polynomial with the given coefficients,
244		* starting from the given initial value.
245		*
246		* @param coefficients the polynomial coefficients array
247		* @param initial the start value to use
248		* @return the point at which the function value is zero
249		* @throws ConvergenceException if the maximum iteration count is exceeded
250		* or the solver detects convergence problems otherwise
251		* @throws FunctionEvaluationException if an error occurs evaluating the
252		* function
253		* @throws IllegalArgumentException if any parameters are invalid
254		*/
255		public Complex solve(Complex coefficients[], Complex initial) throws
256		ConvergenceException, FunctionEvaluationException {
257
258	32	Complex z = initial, oldz, pv, dv, d2v, G, G2, H, delta, denominator;
259
260	32	int n = coefficients.length - 1;
261	32	if (n < 1) {
262	0	throw new IllegalArgumentException
263		("Polynomial degree must be positive: degree=" + n);
264		}
265	32	Complex N = new Complex((double)n, 0.0);
266	32	Complex N1 = new Complex((double)(n-1), 0.0);
267
268	32	int i = 1;
269	32	oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
270	112	while (i <= maximalIterationCount) {
271		// Compute pv (polynomial value), dv (derivative value), and
272		// d2v (second derivative value) simultaneously.
273	112	pv = coefficients[n];
274	112	dv = d2v = new Complex(0.0, 0.0);
275	508	for (int j = n-1; j >= 0; j--) {
276	396	d2v = dv.add(z.multiply(d2v));
277	396	dv = pv.add(z.multiply(dv));
278	396	pv = coefficients[j].add(z.multiply(pv));
279		}
280	112	d2v = d2v.multiply(new Complex(2.0, 0.0));
281
282		// check for convergence
283	112	double tolerance = Math.max(relativeAccuracy * z.abs(),
284		absoluteAccuracy);
285	112	if ((z.subtract(oldz)).abs() <= tolerance) {
286	14	resultComputed = true;
287	14	iterationCount = i;
288	14	return z;
289		}
290	98	if (pv.abs() <= functionValueAccuracy) {
291	18	resultComputed = true;
292	18	iterationCount = i;
293	18	return z;
294		}
295
296		// now pv != 0, calculate the new approximation
297	80	G = dv.divide(pv);
298	80	G2 = G.multiply(G);
299	80	H = G2.subtract(d2v.divide(pv));
300	80	delta = N1.multiply((N.multiply(H)).subtract(G2));
301		// choose a denominator larger in magnitude
302	80	Complex dplus = G.add(ComplexUtils.sqrt(delta));
303	80	Complex dminus = G.subtract(ComplexUtils.sqrt(delta));
304	80	denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
305		// Perturb z if denominator is zero, for instance,
306		// p(x) = x^3 + 1, z = 0.
307	80	if (denominator.equals(new Complex(0.0, 0.0))) {
308	2	z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
309	2	oldz = new Complex(Double.POSITIVE_INFINITY,
310		Double.POSITIVE_INFINITY);
311		} else {
312	78	oldz = z;
313	78	z = z.subtract(N.divide(denominator));
314		}
315	80	i++;
316		}
317	0	throw new ConvergenceException("Maximum number of iterations exceeded.");
318		}
319		}