View Javadoc

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