    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.math3.analysis.solvers;
    018
    019    import org.apache.commons.math3.complex.Complex;
    020    import org.apache.commons.math3.complex.ComplexUtils;
    021    import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
    022    import org.apache.commons.math3.exception.NoBracketingException;
    023    import org.apache.commons.math3.exception.NullArgumentException;
    024    import org.apache.commons.math3.exception.NoDataException;
    025    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    026    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    027    import org.apache.commons.math3.exception.util.LocalizedFormats;
    028    import org.apache.commons.math3.util.FastMath;
    029
    030    /**
    031     * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
    032     * Laguerre's Method</a> for root finding of real coefficient polynomials.
    033     * For reference, see
    034     * <quote>
    035     *  <b>A First Course in Numerical Analysis</b>
    036     *  ISBN 048641454X, chapter 8.
    037     * </quote>
    038     * Laguerre's method is global in the sense that it can start with any initial
    039     * approximation and be able to solve all roots from that point.
    040     * The algorithm requires a bracketing condition.
    041     *
    042     * @version $Id: LaguerreSolver.java 1422195 2012-12-15 06:45:18Z psteitz $
    043     * @since 1.2
    044     */
    045    public class LaguerreSolver extends AbstractPolynomialSolver {
    046        /** Default absolute accuracy. */
    047        private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    048        /** Complex solver. */
    049        private final ComplexSolver complexSolver = new ComplexSolver();
    050
    051        /**
    052         * Construct a solver with default accuracy (1e-6).
    053         */
    054        public LaguerreSolver() {
    055            this(DEFAULT_ABSOLUTE_ACCURACY);
    056        }
    057        /**
    058         * Construct a solver.
    059         *
    060         * @param absoluteAccuracy Absolute accuracy.
    061         */
    062        public LaguerreSolver(double absoluteAccuracy) {
    063            super(absoluteAccuracy);
    064        }
    065        /**
    066         * Construct a solver.
    067         *
    068         * @param relativeAccuracy Relative accuracy.
    069         * @param absoluteAccuracy Absolute accuracy.
    070         */
    071        public LaguerreSolver(double relativeAccuracy,
    072                              double absoluteAccuracy) {
    073            super(relativeAccuracy, absoluteAccuracy);
    074        }
    075        /**
    076         * Construct a solver.
    077         *
    078         * @param relativeAccuracy Relative accuracy.
    079         * @param absoluteAccuracy Absolute accuracy.
    080         * @param functionValueAccuracy Function value accuracy.
    081         */
    082        public LaguerreSolver(double relativeAccuracy,
    083                              double absoluteAccuracy,
    084                              double functionValueAccuracy) {
    085            super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
    086        }
    087
    088        /**
    089         * {@inheritDoc}
    090         */
    091        @Override
    092        public double doSolve()
    093            throws TooManyEvaluationsException,
    094                   NumberIsTooLargeException,
    095                   NoBracketingException {
    096            final double min = getMin();
    097            final double max = getMax();
    098            final double initial = getStartValue();
    099            final double functionValueAccuracy = getFunctionValueAccuracy();
    100
    101            verifySequence(min, initial, max);
    102
    103            // Return the initial guess if it is good enough.
    104            final double yInitial = computeObjectiveValue(initial);
    105            if (FastMath.abs(yInitial) <= functionValueAccuracy) {
    106                return initial;
    107            }
    108
    109            // Return the first endpoint if it is good enough.
    110            final double yMin = computeObjectiveValue(min);
    111            if (FastMath.abs(yMin) <= functionValueAccuracy) {
    112                return min;
    113            }
    114
    115            // Reduce interval if min and initial bracket the root.
    116            if (yInitial * yMin < 0) {
    117                return laguerre(min, initial, yMin, yInitial);
    118            }
    119
    120            // Return the second endpoint if it is good enough.
    121            final double yMax = computeObjectiveValue(max);
    122            if (FastMath.abs(yMax) <= functionValueAccuracy) {
    123                return max;
    124            }
    125
    126            // Reduce interval if initial and max bracket the root.
    127            if (yInitial * yMax < 0) {
    128                return laguerre(initial, max, yInitial, yMax);
    129            }
    130
    131            throw new NoBracketingException(min, max, yMin, yMax);
    132        }
    133
    134        /**
    135         * Find a real root in the given interval.
    136         *
    137         * Despite the bracketing condition, the root returned by
    138         * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
    139         * not be a real zero inside {@code [min, max]}.
    140         * For example, <code>p(x) = x<sup>3</sup> + 1,</code>
    141         * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
    142         * When it occurs, this code calls
    143         * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
    144         * in order to obtain all roots and picks up one real root.
    145         *
    146         * @param lo Lower bound of the search interval.
    147         * @param hi Higher bound of the search interval.
    148         * @param fLo Function value at the lower bound of the search interval.
    149         * @param fHi Function value at the higher bound of the search interval.
    150         * @return the point at which the function value is zero.
    151         * @deprecated This method should not be part of the public API: It will
    152         * be made private in version 4.0.
    153         */
    154        @Deprecated
    155        public double laguerre(double lo, double hi,
    156                               double fLo, double fHi) {
    157            final Complex c[] = ComplexUtils.convertToComplex(getCoefficients());
    158
    159            final Complex initial = new Complex(0.5 * (lo + hi), 0);
    160            final Complex z = complexSolver.solve(c, initial);
    161            if (complexSolver.isRoot(lo, hi, z)) {
    162                return z.getReal();
    163            } else {
    164                double r = Double.NaN;
    165                // Solve all roots and select the one we are seeking.
    166                Complex[] root = complexSolver.solveAll(c, initial);
    167                for (int i = 0; i < root.length; i++) {
    168                    if (complexSolver.isRoot(lo, hi, root[i])) {
    169                        r = root[i].getReal();
    170                        break;
    171                    }
    172                }
    173                return r;
    174            }
    175        }
    176
    177        /**
    178         * Find all complex roots for the polynomial with the given
    179         * coefficients, starting from the given initial value.
    180         * <br/>
    181         * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
    182         *
    183         * @param coefficients Polynomial coefficients.
    184         * @param initial Start value.
    185         * @return the point at which the function value is zero.
    186         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    187         * if the maximum number of evaluations is exceeded.
    188         * @throws NullArgumentException if the {@code coefficients} is
    189         * {@code null}.
    190         * @throws NoDataException if the {@code coefficients} array is empty.
    191         * @since 3.1
    192         */
    193        public Complex[] solveAllComplex(double[] coefficients,
    194                                         double initial)
    195            throws NullArgumentException,
    196                   NoDataException,
    197                   TooManyEvaluationsException {
    198            setup(Integer.MAX_VALUE,
    199                  new PolynomialFunction(coefficients),
    200                  Double.NEGATIVE_INFINITY,
    201                  Double.POSITIVE_INFINITY,
    202                  initial);
    203            return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
    204                                          new Complex(initial, 0d));
    205        }
    206
    207        /**
    208         * Find a complex root for the polynomial with the given coefficients,
    209         * starting from the given initial value.
    210         * <br/>
    211         * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
    212         *
    213         * @param coefficients Polynomial coefficients.
    214         * @param initial Start value.
    215         * @return the point at which the function value is zero.
    216         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    217         * if the maximum number of evaluations is exceeded.
    218         * @throws NullArgumentException if the {@code coefficients} is
    219         * {@code null}.
    220         * @throws NoDataException if the {@code coefficients} array is empty.
    221         * @since 3.1
    222         */
    223        public Complex solveComplex(double[] coefficients,
    224                                    double initial)
    225            throws NullArgumentException,
    226                   NoDataException,
    227                   TooManyEvaluationsException {
    228            setup(Integer.MAX_VALUE,
    229                  new PolynomialFunction(coefficients),
    230                  Double.NEGATIVE_INFINITY,
    231                  Double.POSITIVE_INFINITY,
    232                  initial);
    233            return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
    234                                       new Complex(initial, 0d));
    235        }
    236
    237        /**
    238         * Class for searching all (complex) roots.
    239         */
    240        private class ComplexSolver {
    241            /**
    242             * Check whether the given complex root is actually a real zero
    243             * in the given interval, within the solver tolerance level.
    244             *
    245             * @param min Lower bound for the interval.
    246             * @param max Upper bound for the interval.
    247             * @param z Complex root.
    248             * @return {@code true} if z is a real zero.
    249             */
    250            public boolean isRoot(double min, double max, Complex z) {
    251                if (isSequence(min, z.getReal(), max)) {
    252                    double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
    253                    return (FastMath.abs(z.getImaginary()) <= tolerance) ||
    254                         (z.abs() <= getFunctionValueAccuracy());
    255                }
    256                return false;
    257            }
    258
    259            /**
    260             * Find all complex roots for the polynomial with the given
    261             * coefficients, starting from the given initial value.
    262             *
    263             * @param coefficients Polynomial coefficients.
    264             * @param initial Start value.
    265             * @return the point at which the function value is zero.
    266             * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    267             * if the maximum number of evaluations is exceeded.
    268             * @throws NullArgumentException if the {@code coefficients} is
    269             * {@code null}.
    270             * @throws NoDataException if the {@code coefficients} array is empty.
    271             */
    272            public Complex[] solveAll(Complex coefficients[], Complex initial)
    273                throws NullArgumentException,
    274                       NoDataException,
    275                       TooManyEvaluationsException {
    276                if (coefficients == null) {
    277                    throw new NullArgumentException();
    278                }
    279                final int n = coefficients.length - 1;
    280                if (n == 0) {
    281                    throw new NoDataException(LocalizedFormats.POLYNOMIAL);
    282                }
    283                // Coefficients for deflated polynomial.
    284                final Complex c[] = new Complex[n + 1];
    285                for (int i = 0; i <= n; i++) {
    286                    c[i] = coefficients[i];
    287                }
    288
    289                // Solve individual roots successively.
    290                final Complex root[] = new Complex[n];
    291                for (int i = 0; i < n; i++) {
    292                    final Complex subarray[] = new Complex[n - i + 1];
    293                    System.arraycopy(c, 0, subarray, 0, subarray.length);
    294                    root[i] = solve(subarray, initial);
    295                    // Polynomial deflation using synthetic division.
    296                    Complex newc = c[n - i];
    297                    Complex oldc = null;
    298                    for (int j = n - i - 1; j >= 0; j--) {
    299                        oldc = c[j];
    300                        c[j] = newc;
    301                        newc = oldc.add(newc.multiply(root[i]));
    302                    }
    303                }
    304
    305                return root;
    306            }
    307
    308            /**
    309             * Find a complex root for the polynomial with the given coefficients,
    310             * starting from the given initial value.
    311             *
    312             * @param coefficients Polynomial coefficients.
    313             * @param initial Start value.
    314             * @return the point at which the function value is zero.
    315             * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    316             * if the maximum number of evaluations is exceeded.
    317             * @throws NullArgumentException if the {@code coefficients} is
    318             * {@code null}.
    319             * @throws NoDataException if the {@code coefficients} array is empty.
    320             */
    321            public Complex solve(Complex coefficients[], Complex initial)
    322                throws NullArgumentException,
    323                       NoDataException,
    324                       TooManyEvaluationsException {
    325                if (coefficients == null) {
    326                    throw new NullArgumentException();
    327                }
    328
    329                final int n = coefficients.length - 1;
    330                if (n == 0) {
    331                    throw new NoDataException(LocalizedFormats.POLYNOMIAL);
    332                }
    333
    334                final double absoluteAccuracy = getAbsoluteAccuracy();
    335                final double relativeAccuracy = getRelativeAccuracy();
    336                final double functionValueAccuracy = getFunctionValueAccuracy();
    337
    338                final Complex nC  = new Complex(n, 0);
    339                final Complex n1C = new Complex(n - 1, 0);
    340
    341                Complex z = initial;
    342                Complex oldz = new Complex(Double.POSITIVE_INFINITY,
    343                                           Double.POSITIVE_INFINITY);
    344                while (true) {
    345                    // Compute pv (polynomial value), dv (derivative value), and
    346                    // d2v (second derivative value) simultaneously.
    347                    Complex pv = coefficients[n];
    348                    Complex dv = Complex.ZERO;
    349                    Complex d2v = Complex.ZERO;
    350                    for (int j = n-1; j >= 0; j--) {
    351                        d2v = dv.add(z.multiply(d2v));
    352                        dv = pv.add(z.multiply(dv));
    353                        pv = coefficients[j].add(z.multiply(pv));
    354                    }
    355                    d2v = d2v.multiply(new Complex(2.0, 0.0));
    356
    357                    // Check for convergence.
    358                    final double tolerance = FastMath.max(relativeAccuracy * z.abs(),
    359                                                          absoluteAccuracy);
    360                    if ((z.subtract(oldz)).abs() <= tolerance) {
    361                        return z;
    362                    }
    363                    if (pv.abs() <= functionValueAccuracy) {
    364                        return z;
    365                    }
    366
    367                    // Now pv != 0, calculate the new approximation.
    368                    final Complex G = dv.divide(pv);
    369                    final Complex G2 = G.multiply(G);
    370                    final Complex H = G2.subtract(d2v.divide(pv));
    371                    final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
    372                    // Choose a denominator larger in magnitude.
    373                    final Complex deltaSqrt = delta.sqrt();
    374                    final Complex dplus = G.add(deltaSqrt);
    375                    final Complex dminus = G.subtract(deltaSqrt);
    376                    final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
    377                    // Perturb z if denominator is zero, for instance,
    378                    // p(x) = x^3 + 1, z = 0.
    379                    if (denominator.equals(new Complex(0.0, 0.0))) {
    380                        z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
    381                        oldz = new Complex(Double.POSITIVE_INFINITY,
    382                                           Double.POSITIVE_INFINITY);
    383                    } else {
    384                        oldz = z;
    385                        z = z.subtract(nC.divide(denominator));
    386                    }
    387                    incrementEvaluationCount();
    388                }
    389            }
    390        }
    391    }