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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
037     * @since 1.2
038     */
039    public class LaguerreSolver extends UnivariateRealSolverImpl {
040        /** polynomial function to solve.
041         * @deprecated as of 2.0 the function is not stored anymore in the instance
042         */
043        @Deprecated
044        private PolynomialFunction p;
045    
046        /**
047         * Construct a solver for the given function.
048         *
049         * @param f function to solve
050         * @throws IllegalArgumentException if function is not polynomial
051         * @deprecated as of 2.0 the function to solve is passed as an argument
052         * to the {@link #solve(UnivariateRealFunction, double, double)} or
053         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
054         * method.
055         */
056        @Deprecated
057        public LaguerreSolver(UnivariateRealFunction f) throws
058            IllegalArgumentException {
059            super(f, 100, 1E-6);
060            if (f instanceof PolynomialFunction) {
061                p = (PolynomialFunction) f;
062            } else {
063                throw MathRuntimeException.createIllegalArgumentException("function is not polynomial");
064            }
065        }
066    
067        /**
068         * Construct a solver.
069         */
070        public LaguerreSolver() {
071            super(100, 1E-6);
072        }
073    
074        /**
075         * Returns a copy of the polynomial function.
076         * 
077         * @return a fresh copy of the polynomial function
078         * @deprecated as of 2.0 the function is not stored anymore within the instance.
079         */
080        @Deprecated
081        public PolynomialFunction getPolynomialFunction() {
082            return new PolynomialFunction(p.getCoefficients());
083        }
084    
085        /** {@inheritDoc} */
086        @Deprecated
087        public double solve(final double min, final double max)
088            throws ConvergenceException, FunctionEvaluationException {
089            return solve(p, min, max);
090        }
091    
092        /** {@inheritDoc} */
093        @Deprecated
094        public double solve(final double min, final double max, final double initial)
095            throws ConvergenceException, FunctionEvaluationException {
096            return solve(p, min, max, initial);
097        }
098    
099        /**
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    }