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    import org.apache.commons.math.exception.util.LocalizedFormats;
027    import org.apache.commons.math.util.FastMath;
028    
029    /**
030     * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
031     * Laguerre's Method</a> for root finding of real coefficient polynomials.
032     * For reference, see <b>A First Course in Numerical Analysis</b>,
033     * ISBN 048641454X, chapter 8.
034     * <p>
035     * Laguerre's method is global in the sense that it can start with any initial
036     * approximation and be able to solve all roots from that point.</p>
037     *
038     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
039     * @since 1.2
040     */
041    public class LaguerreSolver extends UnivariateRealSolverImpl {
042    
043        /** polynomial function to solve.
044         * @deprecated as of 2.0 the function is not stored anymore in the instance
045         */
046        @Deprecated
047        private final PolynomialFunction p;
048    
049        /**
050         * Construct a solver for the given function.
051         *
052         * @param f function to solve
053         * @throws IllegalArgumentException if function is not polynomial
054         * @deprecated as of 2.0 the function to solve is passed as an argument
055         * to the {@link #solve(UnivariateRealFunction, double, double)} or
056         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
057         * method.
058         */
059        @Deprecated
060        public LaguerreSolver(UnivariateRealFunction f) throws IllegalArgumentException {
061            super(f, 100, 1E-6);
062            if (f instanceof PolynomialFunction) {
063                p = (PolynomialFunction) f;
064            } else {
065                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL);
066            }
067        }
068    
069        /**
070         * Construct a solver.
071         * @deprecated in 2.2 (to be removed in 3.0)
072         */
073        @Deprecated
074        public LaguerreSolver() {
075            super(100, 1E-6);
076            p = null;
077        }
078    
079        /**
080         * Returns a copy of the polynomial function.
081         *
082         * @return a fresh copy of the polynomial function
083         * @deprecated as of 2.0 the function is not stored anymore within the instance.
084         */
085        @Deprecated
086        public PolynomialFunction getPolynomialFunction() {
087            return new PolynomialFunction(p.getCoefficients());
088        }
089    
090        /** {@inheritDoc} */
091        @Deprecated
092        public double solve(final double min, final double max)
093            throws ConvergenceException, FunctionEvaluationException {
094            return solve(p, min, max);
095        }
096    
097        /** {@inheritDoc} */
098        @Deprecated
099        public double solve(final double min, final double max, final double initial)
100            throws ConvergenceException, FunctionEvaluationException {
101            return solve(p, min, max, initial);
102        }
103    
104        /**
105         * Find a real root in the given interval with initial value.
106         * <p>
107         * Requires bracketing condition.</p>
108         *
109         * @param f function to solve (must be polynomial)
110         * @param min the lower bound for the interval
111         * @param max the upper bound for the interval
112         * @param initial the start value to use
113         * @param maxEval Maximum number of evaluations.
114         * @return the point at which the function value is zero
115         * @throws ConvergenceException if the maximum iteration count is exceeded
116         * or the solver detects convergence problems otherwise
117         * @throws FunctionEvaluationException if an error occurs evaluating the function
118         * @throws IllegalArgumentException if any parameters are invalid
119         */
120        @Override
121        public double solve(int maxEval, final UnivariateRealFunction f,
122                            final double min, final double max, final double initial)
123            throws ConvergenceException, FunctionEvaluationException {
124            setMaximalIterationCount(maxEval);
125            return solve(f, min, max, initial);
126        }
127    
128        /**
129         * Find a real root in the given interval with initial value.
130         * <p>
131         * Requires bracketing condition.</p>
132         *
133         * @param f function to solve (must be polynomial)
134         * @param min the lower bound for the interval
135         * @param max the upper bound for the interval
136         * @param initial the start value to use
137         * @return the point at which the function value is zero
138         * @throws ConvergenceException if the maximum iteration count is exceeded
139         * or the solver detects convergence problems otherwise
140         * @throws FunctionEvaluationException if an error occurs evaluating the function
141         * @throws IllegalArgumentException if any parameters are invalid
142         * @deprecated in 2.2 (to be removed in 3.0).
143         */
144        @Deprecated
145        public double solve(final UnivariateRealFunction f,
146                            final double min, final double max, final double initial)
147            throws ConvergenceException, FunctionEvaluationException {
148    
149            // check for zeros before verifying bracketing
150            if (f.value(min) == 0.0) {
151                return min;
152            }
153            if (f.value(max) == 0.0) {
154                return max;
155            }
156            if (f.value(initial) == 0.0) {
157                return initial;
158            }
159    
160            verifyBracketing(min, max, f);
161            verifySequence(min, initial, max);
162            if (isBracketing(min, initial, f)) {
163                return solve(f, min, initial);
164            } else {
165                return solve(f, initial, max);
166            }
167    
168        }
169    
170        /**
171         * Find a real root in the given interval.
172         * <p>
173         * Despite the bracketing condition, the root returned by solve(Complex[],
174         * Complex) may not be a real zero inside [min, max]. For example,
175         * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
176         * another initial value, or, as we did here, call solveAll() to obtain
177         * all roots and pick up the one that we're looking for.</p>
178         *
179         * @param f the function to solve
180         * @param min the lower bound for the interval
181         * @param max the upper bound for the interval
182         * @param maxEval Maximum number of evaluations.
183         * @return the point at which the function value is zero
184         * @throws ConvergenceException if the maximum iteration count is exceeded
185         * or the solver detects convergence problems otherwise
186         * @throws FunctionEvaluationException if an error occurs evaluating the function
187         * @throws IllegalArgumentException if any parameters are invalid
188         */
189        @Override
190        public double solve(int maxEval, final UnivariateRealFunction f,
191                            final double min, final double max)
192            throws ConvergenceException, FunctionEvaluationException {
193            setMaximalIterationCount(maxEval);
194            return solve(f, min, max);
195        }
196    
197        /**
198         * Find a real root in the given interval.
199         * <p>
200         * Despite the bracketing condition, the root returned by solve(Complex[],
201         * Complex) may not be a real zero inside [min, max]. For example,
202         * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
203         * another initial value, or, as we did here, call solveAll() to obtain
204         * all roots and pick up the one that we're looking for.</p>
205         *
206         * @param f the function to solve
207         * @param min the lower bound for the interval
208         * @param max the upper bound for the interval
209         * @return the point at which the function value is zero
210         * @throws ConvergenceException if the maximum iteration count is exceeded
211         * or the solver detects convergence problems otherwise
212         * @throws FunctionEvaluationException if an error occurs evaluating the function
213         * @throws IllegalArgumentException if any parameters are invalid
214         * @deprecated in 2.2 (to be removed in 3.0).
215         */
216        @Deprecated
217        public double solve(final UnivariateRealFunction f,
218                            final double min, final double max)
219            throws ConvergenceException, FunctionEvaluationException {
220    
221            // check function type
222            if (!(f instanceof PolynomialFunction)) {
223                throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL);
224            }
225    
226            // check for zeros before verifying bracketing
227            if (f.value(min) == 0.0) { return min; }
228            if (f.value(max) == 0.0) { return max; }
229            verifyBracketing(min, max, f);
230    
231            double coefficients[] = ((PolynomialFunction) f).getCoefficients();
232            Complex c[] = new Complex[coefficients.length];
233            for (int i = 0; i < coefficients.length; i++) {
234                c[i] = new Complex(coefficients[i], 0.0);
235            }
236            Complex initial = new Complex(0.5 * (min + max), 0.0);
237            Complex z = solve(c, initial);
238            if (isRootOK(min, max, z)) {
239                setResult(z.getReal(), iterationCount);
240                return result;
241            }
242    
243            // solve all roots and select the one we're seeking
244            Complex[] root = solveAll(c, initial);
245            for (int i = 0; i < root.length; i++) {
246                if (isRootOK(min, max, root[i])) {
247                    setResult(root[i].getReal(), iterationCount);
248                    return result;
249                }
250            }
251    
252            // should never happen
253            throw new ConvergenceException();
254        }
255    
256        /**
257         * Returns true iff the given complex root is actually a real zero
258         * in the given interval, within the solver tolerance level.
259         *
260         * @param min the lower bound for the interval
261         * @param max the upper bound for the interval
262         * @param z the complex root
263         * @return true iff z is the sought-after real zero
264         */
265        protected boolean isRootOK(double min, double max, Complex z) {
266            double tolerance = FastMath.max(relativeAccuracy * z.abs(), absoluteAccuracy);
267            return (isSequence(min, z.getReal(), max)) &&
268                   (FastMath.abs(z.getImaginary()) <= tolerance ||
269                    z.abs() <= functionValueAccuracy);
270        }
271    
272        /**
273         * Find all complex roots for the polynomial with the given coefficients,
274         * starting from the given initial value.
275         *
276         * @param coefficients the polynomial coefficients array
277         * @param initial the start value to use
278         * @return the point at which the function value is zero
279         * @throws ConvergenceException if the maximum iteration count is exceeded
280         * or the solver detects convergence problems otherwise
281         * @throws FunctionEvaluationException if an error occurs evaluating the function
282         * @throws IllegalArgumentException if any parameters are invalid
283         * @deprecated in 2.2.
284         */
285        @Deprecated
286        public Complex[] solveAll(double coefficients[], double initial) throws
287            ConvergenceException, FunctionEvaluationException {
288    
289            Complex c[] = new Complex[coefficients.length];
290            Complex z = new Complex(initial, 0.0);
291            for (int i = 0; i < c.length; i++) {
292                c[i] = new Complex(coefficients[i], 0.0);
293            }
294            return solveAll(c, z);
295        }
296    
297        /**
298         * Find all complex roots 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 function
307         * @throws IllegalArgumentException if any parameters are invalid
308         * @deprecated in 2.2.
309         */
310        @Deprecated
311        public Complex[] solveAll(Complex coefficients[], Complex initial) throws
312            MaxIterationsExceededException, FunctionEvaluationException {
313    
314            int n = coefficients.length - 1;
315            int iterationCount = 0;
316            if (n < 1) {
317                throw MathRuntimeException.createIllegalArgumentException(
318                      LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
319            }
320            Complex c[] = new Complex[n+1];    // coefficients for deflated polynomial
321            for (int i = 0; i <= n; i++) {
322                c[i] = coefficients[i];
323            }
324    
325            // solve individual root successively
326            Complex root[] = new Complex[n];
327            for (int i = 0; i < n; i++) {
328                Complex subarray[] = new Complex[n-i+1];
329                System.arraycopy(c, 0, subarray, 0, subarray.length);
330                root[i] = solve(subarray, initial);
331                // polynomial deflation using synthetic division
332                Complex newc = c[n-i];
333                Complex oldc = null;
334                for (int j = n-i-1; j >= 0; j--) {
335                    oldc = c[j];
336                    c[j] = newc;
337                    newc = oldc.add(newc.multiply(root[i]));
338                }
339                iterationCount += this.iterationCount;
340            }
341    
342            resultComputed = true;
343            this.iterationCount = iterationCount;
344            return root;
345        }
346    
347        /**
348         * Find a complex root for the polynomial with the given coefficients,
349         * starting from the given initial value.
350         *
351         * @param coefficients the polynomial coefficients array
352         * @param initial the start value to use
353         * @return the point at which the function value is zero
354         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
355         * or the solver detects convergence problems otherwise
356         * @throws FunctionEvaluationException if an error occurs evaluating the function
357         * @throws IllegalArgumentException if any parameters are invalid
358         * @deprecated in 2.2.
359         */
360        @Deprecated
361        public Complex solve(Complex coefficients[], Complex initial) throws
362            MaxIterationsExceededException, FunctionEvaluationException {
363    
364            int n = coefficients.length - 1;
365            if (n < 1) {
366                throw MathRuntimeException.createIllegalArgumentException(
367                      LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
368            }
369            Complex N  = new Complex(n,     0.0);
370            Complex N1 = new Complex(n - 1, 0.0);
371    
372            int i = 1;
373            Complex pv = null;
374            Complex dv = null;
375            Complex d2v = null;
376            Complex G = null;
377            Complex G2 = null;
378            Complex H = null;
379            Complex delta = null;
380            Complex denominator = null;
381            Complex z = initial;
382            Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
383            while (i <= maximalIterationCount) {
384                // Compute pv (polynomial value), dv (derivative value), and
385                // d2v (second derivative value) simultaneously.
386                pv = coefficients[n];
387                dv = Complex.ZERO;
388                d2v = Complex.ZERO;
389                for (int j = n-1; j >= 0; j--) {
390                    d2v = dv.add(z.multiply(d2v));
391                    dv = pv.add(z.multiply(dv));
392                    pv = coefficients[j].add(z.multiply(pv));
393                }
394                d2v = d2v.multiply(new Complex(2.0, 0.0));
395    
396                // check for convergence
397                double tolerance = FastMath.max(relativeAccuracy * z.abs(),
398                                            absoluteAccuracy);
399                if ((z.subtract(oldz)).abs() <= tolerance) {
400                    resultComputed = true;
401                    iterationCount = i;
402                    return z;
403                }
404                if (pv.abs() <= functionValueAccuracy) {
405                    resultComputed = true;
406                    iterationCount = i;
407                    return z;
408                }
409    
410                // now pv != 0, calculate the new approximation
411                G = dv.divide(pv);
412                G2 = G.multiply(G);
413                H = G2.subtract(d2v.divide(pv));
414                delta = N1.multiply((N.multiply(H)).subtract(G2));
415                // choose a denominator larger in magnitude
416                Complex deltaSqrt = delta.sqrt();
417                Complex dplus = G.add(deltaSqrt);
418                Complex dminus = G.subtract(deltaSqrt);
419                denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
420                // Perturb z if denominator is zero, for instance,
421                // p(x) = x^3 + 1, z = 0.
422                if (denominator.equals(new Complex(0.0, 0.0))) {
423                    z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
424                    oldz = new Complex(Double.POSITIVE_INFINITY,
425                                       Double.POSITIVE_INFINITY);
426                } else {
427                    oldz = z;
428                    z = z.subtract(N.divide(denominator));
429                }
430                i++;
431            }
432            throw new MaxIterationsExceededException(maximalIterationCount);
433        }
434    }