001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.analysis.solvers;
018
019 import org.apache.commons.math.ConvergenceException;
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.MathRuntimeException;
022 import org.apache.commons.math.MaxIterationsExceededException;
023 import org.apache.commons.math.analysis.UnivariateRealFunction;
024 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
025 import org.apache.commons.math.complex.Complex;
026
027 /**
028 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
029 * Laguerre's Method</a> for root finding of real coefficient polynomials.
030 * For reference, see <b>A First Course in Numerical Analysis</b>,
031 * ISBN 048641454X, chapter 8.
032 * <p>
033 * Laguerre's method is global in the sense that it can start with any initial
034 * approximation and be able to solve all roots from that point.</p>
035 *
036 * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $
037 * @since 1.2
038 */
039 public class LaguerreSolver extends UnivariateRealSolverImpl {
040
041 /** Message for non-polynomial function. */
042 private static final String NON_POLYNOMIAL_FUNCTION_MESSAGE =
043 "function is not polynomial";
044
045 /** Message for non-positive degree. */
046 private static final String NON_POSITIVE_DEGREE_MESSAGE =
047 "polynomial degree must be positive: degree={0}";
048
049 /** polynomial function to solve.
050 * @deprecated as of 2.0 the function is not stored anymore in the instance
051 */
052 @Deprecated
053 private final PolynomialFunction p;
054
055 /**
056 * Construct a solver for the given function.
057 *
058 * @param f function to solve
059 * @throws IllegalArgumentException if function is not polynomial
060 * @deprecated as of 2.0 the function to solve is passed as an argument
061 * to the {@link #solve(UnivariateRealFunction, double, double)} or
062 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
063 * method.
064 */
065 @Deprecated
066 public LaguerreSolver(UnivariateRealFunction f) throws
067 IllegalArgumentException {
068 super(f, 100, 1E-6);
069 if (f instanceof PolynomialFunction) {
070 p = (PolynomialFunction) f;
071 } else {
072 throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE);
073 }
074 }
075
076 /**
077 * Construct a solver.
078 */
079 public LaguerreSolver() {
080 super(100, 1E-6);
081 p = null;
082 }
083
084 /**
085 * Returns a copy of the polynomial function.
086 *
087 * @return a fresh copy of the polynomial function
088 * @deprecated as of 2.0 the function is not stored anymore within the instance.
089 */
090 @Deprecated
091 public PolynomialFunction getPolynomialFunction() {
092 return new PolynomialFunction(p.getCoefficients());
093 }
094
095 /** {@inheritDoc} */
096 @Deprecated
097 public double solve(final double min, final double max)
098 throws ConvergenceException, FunctionEvaluationException {
099 return solve(p, min, max);
100 }
101
102 /** {@inheritDoc} */
103 @Deprecated
104 public double solve(final double min, final double max, final double initial)
105 throws ConvergenceException, FunctionEvaluationException {
106 return solve(p, min, max, initial);
107 }
108
109 /**
110 * Find a real root in the given interval with initial value.
111 * <p>
112 * Requires bracketing condition.</p>
113 *
114 * @param f function to solve (must be polynomial)
115 * @param min the lower bound for the interval
116 * @param max the upper bound for the interval
117 * @param initial the start value to use
118 * @return the point at which the function value is zero
119 * @throws ConvergenceException if the maximum iteration count is exceeded
120 * or the solver detects convergence problems otherwise
121 * @throws FunctionEvaluationException if an error occurs evaluating the
122 * function
123 * @throws IllegalArgumentException if any parameters are invalid
124 */
125 public double solve(final UnivariateRealFunction f,
126 final double min, final double max, final double initial)
127 throws ConvergenceException, FunctionEvaluationException {
128
129 // check for zeros before verifying bracketing
130 if (f.value(min) == 0.0) {
131 return min;
132 }
133 if (f.value(max) == 0.0) {
134 return max;
135 }
136 if (f.value(initial) == 0.0) {
137 return initial;
138 }
139
140 verifyBracketing(min, max, f);
141 verifySequence(min, initial, max);
142 if (isBracketing(min, initial, f)) {
143 return solve(f, min, initial);
144 } else {
145 return solve(f, initial, max);
146 }
147
148 }
149
150 /**
151 * Find a real root in the given interval.
152 * <p>
153 * Despite the bracketing condition, the root returned by solve(Complex[],
154 * Complex) may not be a real zero inside [min, max]. For example,
155 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
156 * another initial value, or, as we did here, call solveAll() to obtain
157 * all roots and pick up the one that we're looking for.</p>
158 *
159 * @param f the function to solve
160 * @param min the lower bound for the interval
161 * @param max the upper bound for the interval
162 * @return the point at which the function value is zero
163 * @throws ConvergenceException if the maximum iteration count is exceeded
164 * or the solver detects convergence problems otherwise
165 * @throws FunctionEvaluationException if an error occurs evaluating the
166 * function
167 * @throws IllegalArgumentException if any parameters are invalid
168 */
169 public double solve(final UnivariateRealFunction f,
170 final double min, final double max)
171 throws ConvergenceException, FunctionEvaluationException {
172
173 // check function type
174 if (!(f instanceof PolynomialFunction)) {
175 throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE);
176 }
177
178 // check for zeros before verifying bracketing
179 if (f.value(min) == 0.0) { return min; }
180 if (f.value(max) == 0.0) { return max; }
181 verifyBracketing(min, max, f);
182
183 double coefficients[] = ((PolynomialFunction) f).getCoefficients();
184 Complex c[] = new Complex[coefficients.length];
185 for (int i = 0; i < coefficients.length; i++) {
186 c[i] = new Complex(coefficients[i], 0.0);
187 }
188 Complex initial = new Complex(0.5 * (min + max), 0.0);
189 Complex z = solve(c, initial);
190 if (isRootOK(min, max, z)) {
191 setResult(z.getReal(), iterationCount);
192 return result;
193 }
194
195 // solve all roots and select the one we're seeking
196 Complex[] root = solveAll(c, initial);
197 for (int i = 0; i < root.length; i++) {
198 if (isRootOK(min, max, root[i])) {
199 setResult(root[i].getReal(), iterationCount);
200 return result;
201 }
202 }
203
204 // should never happen
205 throw new ConvergenceException();
206 }
207
208 /**
209 * Returns true iff the given complex root is actually a real zero
210 * in the given interval, within the solver tolerance level.
211 *
212 * @param min the lower bound for the interval
213 * @param max the upper bound for the interval
214 * @param z the complex root
215 * @return true iff z is the sought-after real zero
216 */
217 protected boolean isRootOK(double min, double max, Complex z) {
218 double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
219 return (isSequence(min, z.getReal(), max)) &&
220 (Math.abs(z.getImaginary()) <= tolerance ||
221 z.abs() <= functionValueAccuracy);
222 }
223
224 /**
225 * Find all complex roots for the polynomial with the given coefficients,
226 * starting from the given initial value.
227 *
228 * @param coefficients the polynomial coefficients array
229 * @param initial the start value to use
230 * @return the point at which the function value is zero
231 * @throws ConvergenceException if the maximum iteration count is exceeded
232 * or the solver detects convergence problems otherwise
233 * @throws FunctionEvaluationException if an error occurs evaluating the
234 * function
235 * @throws IllegalArgumentException if any parameters are invalid
236 */
237 public Complex[] solveAll(double coefficients[], double initial) throws
238 ConvergenceException, FunctionEvaluationException {
239
240 Complex c[] = new Complex[coefficients.length];
241 Complex z = new Complex(initial, 0.0);
242 for (int i = 0; i < c.length; i++) {
243 c[i] = new Complex(coefficients[i], 0.0);
244 }
245 return solveAll(c, z);
246 }
247
248 /**
249 * Find all complex roots for the polynomial with the given coefficients,
250 * starting from the given initial value.
251 *
252 * @param coefficients the polynomial coefficients array
253 * @param initial the start value to use
254 * @return the point at which the function value is zero
255 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
256 * or the solver detects convergence problems otherwise
257 * @throws FunctionEvaluationException if an error occurs evaluating the
258 * function
259 * @throws IllegalArgumentException if any parameters are invalid
260 */
261 public Complex[] solveAll(Complex coefficients[], Complex initial) throws
262 MaxIterationsExceededException, FunctionEvaluationException {
263
264 int n = coefficients.length - 1;
265 int iterationCount = 0;
266 if (n < 1) {
267 throw MathRuntimeException.createIllegalArgumentException(
268 NON_POSITIVE_DEGREE_MESSAGE, n);
269 }
270 Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial
271 for (int i = 0; i <= n; i++) {
272 c[i] = coefficients[i];
273 }
274
275 // solve individual root successively
276 Complex root[] = new Complex[n];
277 for (int i = 0; i < n; i++) {
278 Complex subarray[] = new Complex[n-i+1];
279 System.arraycopy(c, 0, subarray, 0, subarray.length);
280 root[i] = solve(subarray, initial);
281 // polynomial deflation using synthetic division
282 Complex newc = c[n-i];
283 Complex oldc = null;
284 for (int j = n-i-1; j >= 0; j--) {
285 oldc = c[j];
286 c[j] = newc;
287 newc = oldc.add(newc.multiply(root[i]));
288 }
289 iterationCount += this.iterationCount;
290 }
291
292 resultComputed = true;
293 this.iterationCount = iterationCount;
294 return root;
295 }
296
297 /**
298 * Find a complex root for the polynomial with the given coefficients,
299 * starting from the given initial value.
300 *
301 * @param coefficients the polynomial coefficients array
302 * @param initial the start value to use
303 * @return the point at which the function value is zero
304 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
305 * or the solver detects convergence problems otherwise
306 * @throws FunctionEvaluationException if an error occurs evaluating the
307 * function
308 * @throws IllegalArgumentException if any parameters are invalid
309 */
310 public Complex solve(Complex coefficients[], Complex initial) throws
311 MaxIterationsExceededException, FunctionEvaluationException {
312
313 int n = coefficients.length - 1;
314 if (n < 1) {
315 throw MathRuntimeException.createIllegalArgumentException(
316 NON_POSITIVE_DEGREE_MESSAGE, n);
317 }
318 Complex N = new Complex(n, 0.0);
319 Complex N1 = new Complex(n - 1, 0.0);
320
321 int i = 1;
322 Complex pv = null;
323 Complex dv = null;
324 Complex d2v = null;
325 Complex G = null;
326 Complex G2 = null;
327 Complex H = null;
328 Complex delta = null;
329 Complex denominator = null;
330 Complex z = initial;
331 Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
332 while (i <= maximalIterationCount) {
333 // Compute pv (polynomial value), dv (derivative value), and
334 // d2v (second derivative value) simultaneously.
335 pv = coefficients[n];
336 dv = Complex.ZERO;
337 d2v = Complex.ZERO;
338 for (int j = n-1; j >= 0; j--) {
339 d2v = dv.add(z.multiply(d2v));
340 dv = pv.add(z.multiply(dv));
341 pv = coefficients[j].add(z.multiply(pv));
342 }
343 d2v = d2v.multiply(new Complex(2.0, 0.0));
344
345 // check for convergence
346 double tolerance = Math.max(relativeAccuracy * z.abs(),
347 absoluteAccuracy);
348 if ((z.subtract(oldz)).abs() <= tolerance) {
349 resultComputed = true;
350 iterationCount = i;
351 return z;
352 }
353 if (pv.abs() <= functionValueAccuracy) {
354 resultComputed = true;
355 iterationCount = i;
356 return z;
357 }
358
359 // now pv != 0, calculate the new approximation
360 G = dv.divide(pv);
361 G2 = G.multiply(G);
362 H = G2.subtract(d2v.divide(pv));
363 delta = N1.multiply((N.multiply(H)).subtract(G2));
364 // choose a denominator larger in magnitude
365 Complex deltaSqrt = delta.sqrt();
366 Complex dplus = G.add(deltaSqrt);
367 Complex dminus = G.subtract(deltaSqrt);
368 denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
369 // Perturb z if denominator is zero, for instance,
370 // p(x) = x^3 + 1, z = 0.
371 if (denominator.equals(new Complex(0.0, 0.0))) {
372 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
373 oldz = new Complex(Double.POSITIVE_INFINITY,
374 Double.POSITIVE_INFINITY);
375 } else {
376 oldz = z;
377 z = z.subtract(N.divide(denominator));
378 }
379 i++;
380 }
381 throw new MaxIterationsExceededException(maximalIterationCount);
382 }
383 }