    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    }