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 }