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    }