View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.solvers;
23  
24  import org.hipparchus.analysis.polynomials.PolynomialFunction;
25  import org.hipparchus.complex.Complex;
26  import org.hipparchus.complex.ComplexUtils;
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.exception.MathIllegalStateException;
30  import org.hipparchus.exception.NullArgumentException;
31  import org.hipparchus.util.FastMath;
32  
33  /**
34   * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
35   * Laguerre's Method</a> for root finding of real coefficient polynomials.
36   * For reference, see
37   * <blockquote>
38   *  <b>A First Course in Numerical Analysis</b>,
39   *  ISBN 048641454X, chapter 8.
40   * </blockquote>
41   * Laguerre's method is global in the sense that it can start with any initial
42   * approximation and be able to solve all roots from that point.
43   * The algorithm requires a bracketing condition.
44   *
45   */
46  public class LaguerreSolver extends AbstractPolynomialSolver {
47      /** Default absolute accuracy. */
48      private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
49      /** Complex solver. */
50      private final ComplexSolver complexSolver = new ComplexSolver();
51  
52      /**
53       * Construct a solver with default accuracy (1e-6).
54       */
55      public LaguerreSolver() {
56          this(DEFAULT_ABSOLUTE_ACCURACY);
57      }
58      /**
59       * Construct a solver.
60       *
61       * @param absoluteAccuracy Absolute accuracy.
62       */
63      public LaguerreSolver(double absoluteAccuracy) {
64          super(absoluteAccuracy);
65      }
66      /**
67       * Construct a solver.
68       *
69       * @param relativeAccuracy Relative accuracy.
70       * @param absoluteAccuracy Absolute accuracy.
71       */
72      public LaguerreSolver(double relativeAccuracy,
73                            double absoluteAccuracy) {
74          super(relativeAccuracy, absoluteAccuracy);
75      }
76      /**
77       * Construct a solver.
78       *
79       * @param relativeAccuracy Relative accuracy.
80       * @param absoluteAccuracy Absolute accuracy.
81       * @param functionValueAccuracy Function value accuracy.
82       */
83      public LaguerreSolver(double relativeAccuracy,
84                            double absoluteAccuracy,
85                            double functionValueAccuracy) {
86          super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
87      }
88  
89      /**
90       * {@inheritDoc}
91       */
92      @Override
93      public double doSolve()
94          throws MathIllegalArgumentException, MathIllegalStateException {
95          final double min = getMin();
96          final double max = getMax();
97          final double initial = getStartValue();
98          final double functionValueAccuracy = getFunctionValueAccuracy();
99  
100         verifySequence(min, initial, max);
101 
102         // Return the initial guess if it is good enough.
103         final double yInitial = computeObjectiveValue(initial);
104         if (FastMath.abs(yInitial) <= functionValueAccuracy) {
105             return initial;
106         }
107 
108         // Return the first endpoint if it is good enough.
109         final double yMin = computeObjectiveValue(min);
110         if (FastMath.abs(yMin) <= functionValueAccuracy) {
111             return min;
112         }
113 
114         // Reduce interval if min and initial bracket the root.
115         if (yInitial * yMin < 0) {
116             return laguerre(min, initial);
117         }
118 
119         // Return the second endpoint if it is good enough.
120         final double yMax = computeObjectiveValue(max);
121         if (FastMath.abs(yMax) <= functionValueAccuracy) {
122             return max;
123         }
124 
125         // Reduce interval if initial and max bracket the root.
126         if (yInitial * yMax < 0) {
127             return laguerre(initial, max);
128         }
129 
130         throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
131                                                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      * @return the point at which the function value is zero.
149      */
150     private double laguerre(double lo, double hi) {
151         final Complex c[] = ComplexUtils.convertToComplex(getCoefficients());
152 
153         final Complex initial = new Complex(0.5 * (lo + hi), 0);
154         final Complex z = complexSolver.solve(c, initial);
155         if (complexSolver.isRoot(lo, hi, z)) {
156             return z.getReal();
157         } else {
158             double r = Double.NaN;
159             // Solve all roots and select the one we are seeking.
160             Complex[] root = complexSolver.solveAll(c, initial);
161             for (int i = 0; i < root.length; i++) {
162                 if (complexSolver.isRoot(lo, hi, root[i])) {
163                     r = root[i].getReal();
164                     break;
165                 }
166             }
167             return r;
168         }
169     }
170 
171     /**
172      * Find all complex roots for the polynomial with the given
173      * coefficients, starting from the given initial value.
174      * <p>
175      * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
176      *
177      * @param coefficients Polynomial coefficients.
178      * @param initial Start value.
179      * @return the point at which the function value is zero.
180      * @throws org.hipparchus.exception.MathIllegalStateException
181      * if the maximum number of evaluations is exceeded.
182      * @throws NullArgumentException if the {@code coefficients} is
183      * {@code null}.
184      * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
185      */
186     public Complex[] solveAllComplex(double[] coefficients,
187                                      double initial)
188         throws MathIllegalArgumentException, NullArgumentException,
189                MathIllegalStateException {
190         setup(Integer.MAX_VALUE,
191               new PolynomialFunction(coefficients),
192               Double.NEGATIVE_INFINITY,
193               Double.POSITIVE_INFINITY,
194               initial);
195         return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
196                                       new Complex(initial, 0d));
197     }
198 
199     /**
200      * Find a complex root for the polynomial with the given coefficients,
201      * starting from the given initial value.
202      * <p>
203      * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
204      *
205      * @param coefficients Polynomial coefficients.
206      * @param initial Start value.
207      * @return the point at which the function value is zero.
208      * @throws org.hipparchus.exception.MathIllegalStateException
209      * if the maximum number of evaluations is exceeded.
210      * @throws NullArgumentException if the {@code coefficients} is
211      * {@code null}.
212      * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
213      */
214     public Complex solveComplex(double[] coefficients,
215                                 double initial)
216         throws MathIllegalArgumentException, NullArgumentException,
217                MathIllegalStateException {
218         setup(Integer.MAX_VALUE,
219               new PolynomialFunction(coefficients),
220               Double.NEGATIVE_INFINITY,
221               Double.POSITIVE_INFINITY,
222               initial);
223         return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
224                                    new Complex(initial, 0d));
225     }
226 
227     /**
228      * Class for searching all (complex) roots.
229      */
230     private class ComplexSolver {
231         /**
232          * Check whether the given complex root is actually a real zero
233          * in the given interval, within the solver tolerance level.
234          *
235          * @param min Lower bound for the interval.
236          * @param max Upper bound for the interval.
237          * @param z Complex root.
238          * @return {@code true} if z is a real zero.
239          */
240         public boolean isRoot(double min, double max, Complex z) {
241             if (isSequence(min, z.getReal(), max)) {
242                 final double zAbs = z.norm();
243                 double tolerance = FastMath.max(getRelativeAccuracy() * zAbs, getAbsoluteAccuracy());
244                 return (FastMath.abs(z.getImaginary()) <= tolerance) ||
245                      (zAbs <= getFunctionValueAccuracy());
246             }
247             return false;
248         }
249 
250         /**
251          * Find all complex roots for the polynomial with the given
252          * coefficients, starting from the given initial value.
253          *
254          * @param coefficients Polynomial coefficients.
255          * @param initial Start value.
256          * @return the point at which the function value is zero.
257          * @throws org.hipparchus.exception.MathIllegalStateException
258          * if the maximum number of evaluations is exceeded.
259          * @throws NullArgumentException if the {@code coefficients} is
260          * {@code null}.
261          * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
262          */
263         public Complex[] solveAll(Complex coefficients[], Complex initial)
264             throws MathIllegalArgumentException, NullArgumentException,
265                    MathIllegalStateException {
266             if (coefficients == null) {
267                 throw new NullArgumentException();
268             }
269             final int n = coefficients.length - 1;
270             if (n == 0) {
271                 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
272             }
273             // Coefficients for deflated polynomial.
274             final Complex c[] = new Complex[n + 1];
275             for (int i = 0; i <= n; i++) {
276                 c[i] = coefficients[i];
277             }
278 
279             // Solve individual roots successively.
280             final Complex root[] = new Complex[n];
281             for (int i = 0; i < n; i++) {
282                 final Complex subarray[] = new Complex[n - i + 1];
283                 System.arraycopy(c, 0, subarray, 0, subarray.length);
284                 root[i] = solve(subarray, initial);
285                 // Polynomial deflation using synthetic division.
286                 Complex newc = c[n - i];
287                 Complex oldc;
288                 for (int j = n - i - 1; j >= 0; j--) {
289                     oldc = c[j];
290                     c[j] = newc;
291                     newc = oldc.add(newc.multiply(root[i]));
292                 }
293             }
294 
295             return root;
296         }
297 
298         /**
299          * Find a complex root for the polynomial with the given coefficients,
300          * starting from the given initial value.
301          *
302          * @param coefficients Polynomial coefficients.
303          * @param initial Start value.
304          * @return the point at which the function value is zero.
305          * @throws org.hipparchus.exception.MathIllegalStateException
306          * if the maximum number of evaluations is exceeded.
307          * @throws NullArgumentException if the {@code coefficients} is
308          * {@code null}.
309          * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
310          */
311         public Complex solve(Complex coefficients[], Complex initial)
312             throws MathIllegalArgumentException, NullArgumentException,
313                    MathIllegalStateException {
314             if (coefficients == null) {
315                 throw new NullArgumentException();
316             }
317 
318             final int n = coefficients.length - 1;
319             if (n == 0) {
320                 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
321             }
322 
323             final double absoluteAccuracy = getAbsoluteAccuracy();
324             final double relativeAccuracy = getRelativeAccuracy();
325             final double functionValueAccuracy = getFunctionValueAccuracy();
326 
327             final Complex nC  = new Complex(n, 0);
328             final Complex n1C = new Complex(n - 1, 0);
329 
330             Complex z = initial;
331             Complex oldz = new Complex(Double.POSITIVE_INFINITY,
332                                        Double.POSITIVE_INFINITY);
333             while (true) {
334                 // Compute pv (polynomial value), dv (derivative value), and
335                 // d2v (second derivative value) simultaneously.
336                 Complex pv = coefficients[n];
337                 Complex dv = Complex.ZERO;
338                 Complex d2v = Complex.ZERO;
339                 for (int j = n-1; j >= 0; j--) {
340                     d2v = dv.add(z.multiply(d2v));
341                     dv = pv.add(z.multiply(dv));
342                     pv = coefficients[j].add(z.multiply(pv));
343                 }
344                 d2v = d2v.multiply(new Complex(2.0, 0.0));
345 
346                 // Check for convergence.
347                 final double tolerance = FastMath.max(relativeAccuracy * z.norm(),
348                                                       absoluteAccuracy);
349                 if ((z.subtract(oldz)).norm() <= tolerance) {
350                     return z;
351                 }
352                 if (pv.norm() <= functionValueAccuracy) {
353                     return z;
354                 }
355 
356                 // Now pv != 0, calculate the new approximation.
357                 final Complex G = dv.divide(pv);
358                 final Complex G2 = G.multiply(G);
359                 final Complex H = G2.subtract(d2v.divide(pv));
360                 final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
361                 // Choose a denominator larger in magnitude.
362                 final Complex deltaSqrt = delta.sqrt();
363                 final Complex dplus = G.add(deltaSqrt);
364                 final Complex dminus = G.subtract(deltaSqrt);
365                 final Complex denominator = dplus.norm() > dminus.norm() ? dplus : dminus;
366                 // Perturb z if denominator is zero, for instance,
367                 // p(x) = x^3 + 1, z = 0.
368                 if (denominator.isZero()) {
369                     z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
370                     oldz = new Complex(Double.POSITIVE_INFINITY,
371                                        Double.POSITIVE_INFINITY);
372                 } else {
373                     oldz = z;
374                     z = z.subtract(nC.divide(denominator));
375                 }
376                 incrementEvaluationCount();
377             }
378         }
379     }
380 }