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    }