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 (Complex complex : root) {
162 if (complexSolver.isRoot(lo, hi, complex)) {
163 r = complex.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 return solveAllComplex(coefficients, 100_000, initial);
191 }
192
193 /**
194 * Find all complex roots for the polynomial with the given
195 * coefficients, starting from the given initial value.
196 * <p>
197 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
198 *
199 * @param coefficients Polynomial coefficients.
200 * @param maxEval Maximum number of evaluations.
201 * @param initial Start value.
202 * @return the point at which the function value is zero.
203 * @throws org.hipparchus.exception.MathIllegalStateException
204 * if the maximum number of evaluations is exceeded.
205 * @throws NullArgumentException if the {@code coefficients} is
206 * {@code null}.
207 * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
208 */
209 public Complex[] solveAllComplex(double[] coefficients,
210 int maxEval,
211 double initial)
212 throws MathIllegalArgumentException, NullArgumentException,
213 MathIllegalStateException {
214 setup(maxEval,
215 new PolynomialFunction(coefficients),
216 Double.NEGATIVE_INFINITY,
217 Double.POSITIVE_INFINITY,
218 initial);
219 return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
220 new Complex(initial, 0d));
221 }
222
223 /**
224 * Find a complex root for the polynomial with the given coefficients,
225 * starting from the given initial value.
226 * <p>
227 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
228 *
229 * @param coefficients Polynomial coefficients.
230 * @param initial Start value.
231 * @return the point at which the function value is zero.
232 * @throws org.hipparchus.exception.MathIllegalStateException
233 * if the maximum number of evaluations is exceeded.
234 * @throws NullArgumentException if the {@code coefficients} is
235 * {@code null}.
236 * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
237 */
238 public Complex solveComplex(double[] coefficients,
239 double initial)
240 throws MathIllegalArgumentException, NullArgumentException,
241 MathIllegalStateException {
242 setup(Integer.MAX_VALUE,
243 new PolynomialFunction(coefficients),
244 Double.NEGATIVE_INFINITY,
245 Double.POSITIVE_INFINITY,
246 initial);
247 return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
248 new Complex(initial, 0d));
249 }
250
251 /**
252 * Class for searching all (complex) roots.
253 */
254 private class ComplexSolver {
255 /**
256 * Check whether the given complex root is actually a real zero
257 * in the given interval, within the solver tolerance level.
258 *
259 * @param min Lower bound for the interval.
260 * @param max Upper bound for the interval.
261 * @param z Complex root.
262 * @return {@code true} if z is a real zero.
263 */
264 public boolean isRoot(double min, double max, Complex z) {
265 if (isSequence(min, z.getReal(), max)) {
266 final double zAbs = z.norm();
267 double tolerance = FastMath.max(getRelativeAccuracy() * zAbs, getAbsoluteAccuracy());
268 return (FastMath.abs(z.getImaginary()) <= tolerance) ||
269 (zAbs <= getFunctionValueAccuracy());
270 }
271 return false;
272 }
273
274 /**
275 * Find all complex roots for the polynomial with the given
276 * coefficients, starting from the given initial value.
277 *
278 * @param coefficients Polynomial coefficients.
279 * @param initial Start value.
280 * @return the point at which the function value is zero.
281 * @throws org.hipparchus.exception.MathIllegalStateException
282 * if the maximum number of evaluations is exceeded.
283 * @throws NullArgumentException if the {@code coefficients} is
284 * {@code null}.
285 * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
286 */
287 public Complex[] solveAll(Complex[] coefficients, Complex initial)
288 throws MathIllegalArgumentException, NullArgumentException,
289 MathIllegalStateException {
290 if (coefficients == null) {
291 throw new NullArgumentException();
292 }
293 final int n = coefficients.length - 1;
294 if (n == 0) {
295 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
296 }
297 // Coefficients for deflated polynomial.
298 final Complex[] c = new Complex[n + 1];
299 System.arraycopy(coefficients, 0, c, 0, n + 1);
300
301 // Solve individual roots successively.
302 final Complex[] root = new Complex[n];
303 for (int i = 0; i < n; i++) {
304 final Complex[] subarray = new Complex[n - i + 1];
305 System.arraycopy(c, 0, subarray, 0, subarray.length);
306 root[i] = solve(subarray, initial);
307 // Polynomial deflation using synthetic division.
308 Complex newc = c[n - i];
309 Complex oldc;
310 for (int j = n - i - 1; j >= 0; j--) {
311 oldc = c[j];
312 c[j] = newc;
313 newc = oldc.add(newc.multiply(root[i]));
314 }
315 }
316
317 return root;
318 }
319
320 /**
321 * Find a complex root for the polynomial with the given coefficients,
322 * starting from the given initial value.
323 *
324 * @param coefficients Polynomial coefficients.
325 * @param initial Start value.
326 * @return the point at which the function value is zero.
327 * @throws org.hipparchus.exception.MathIllegalStateException
328 * if the maximum number of evaluations is exceeded.
329 * @throws NullArgumentException if the {@code coefficients} is
330 * {@code null}.
331 * @throws MathIllegalArgumentException if the {@code coefficients} array is empty.
332 */
333 public Complex solve(Complex[] coefficients, Complex initial)
334 throws MathIllegalArgumentException, NullArgumentException,
335 MathIllegalStateException {
336 if (coefficients == null) {
337 throw new NullArgumentException();
338 }
339
340 final int n = coefficients.length - 1;
341 if (n == 0) {
342 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
343 }
344
345 final double absoluteAccuracy = getAbsoluteAccuracy();
346 final double relativeAccuracy = getRelativeAccuracy();
347 final double functionValueAccuracy = getFunctionValueAccuracy();
348
349 final Complex nC = new Complex(n, 0);
350 final Complex n1C = new Complex(n - 1, 0);
351
352 Complex z = initial;
353 Complex oldz = new Complex(Double.POSITIVE_INFINITY,
354 Double.POSITIVE_INFINITY);
355 while (true) {
356 // Compute pv (polynomial value), dv (derivative value), and
357 // d2v (second derivative value) simultaneously.
358 Complex pv = coefficients[n];
359 Complex dv = Complex.ZERO;
360 Complex d2v = Complex.ZERO;
361 for (int j = n-1; j >= 0; j--) {
362 d2v = dv.add(z.multiply(d2v));
363 dv = pv.add(z.multiply(dv));
364 pv = coefficients[j].add(z.multiply(pv));
365 }
366 d2v = d2v.multiply(new Complex(2.0, 0.0));
367
368 // Check for convergence.
369 final double tolerance = FastMath.max(relativeAccuracy * z.norm(),
370 absoluteAccuracy);
371 if ((z.subtract(oldz)).norm() <= tolerance) {
372 return z;
373 }
374 if (pv.norm() <= functionValueAccuracy) {
375 return z;
376 }
377
378 // Now pv != 0, calculate the new approximation.
379 final Complex G = dv.divide(pv);
380 final Complex G2 = G.multiply(G);
381 final Complex H = G2.subtract(d2v.divide(pv));
382 final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
383 // Choose a denominator larger in magnitude.
384 final Complex deltaSqrt = delta.sqrt();
385 final Complex dplus = G.add(deltaSqrt);
386 final Complex dminus = G.subtract(deltaSqrt);
387 final Complex denominator = dplus.norm() > dminus.norm() ? dplus : dminus;
388 // Perturb z if denominator is zero, for instance,
389 // p(x) = x^3 + 1, z = 0.
390 if (denominator.isZero()) {
391 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
392 oldz = new Complex(Double.POSITIVE_INFINITY,
393 Double.POSITIVE_INFINITY);
394 } else {
395 oldz = z;
396 z = z.subtract(nC.divide(denominator));
397 }
398 incrementEvaluationCount();
399 }
400 }
401 }
402 }