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 */ 017package org.apache.commons.math3.analysis.solvers; 018 019import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; 020import org.apache.commons.math3.complex.Complex; 021import org.apache.commons.math3.complex.ComplexUtils; 022import org.apache.commons.math3.exception.NoBracketingException; 023import org.apache.commons.math3.exception.NoDataException; 024import org.apache.commons.math3.exception.NullArgumentException; 025import org.apache.commons.math3.exception.NumberIsTooLargeException; 026import org.apache.commons.math3.exception.TooManyEvaluationsException; 027import org.apache.commons.math3.exception.util.LocalizedFormats; 028import 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 * <blockquote> 035 * <b>A First Course in Numerical Analysis</b>, 036 * ISBN 048641454X, chapter 8. 037 * </blockquote> 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 * @since 1.2 043 */ 044public class LaguerreSolver extends AbstractPolynomialSolver { 045 /** Default absolute accuracy. */ 046 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; 047 /** Complex solver. */ 048 private final ComplexSolver complexSolver = new ComplexSolver(); 049 050 /** 051 * Construct a solver with default accuracy (1e-6). 052 */ 053 public LaguerreSolver() { 054 this(DEFAULT_ABSOLUTE_ACCURACY); 055 } 056 /** 057 * Construct a solver. 058 * 059 * @param absoluteAccuracy Absolute accuracy. 060 */ 061 public LaguerreSolver(double absoluteAccuracy) { 062 super(absoluteAccuracy); 063 } 064 /** 065 * Construct a solver. 066 * 067 * @param relativeAccuracy Relative accuracy. 068 * @param absoluteAccuracy Absolute accuracy. 069 */ 070 public LaguerreSolver(double relativeAccuracy, 071 double absoluteAccuracy) { 072 super(relativeAccuracy, absoluteAccuracy); 073 } 074 /** 075 * Construct a solver. 076 * 077 * @param relativeAccuracy Relative accuracy. 078 * @param absoluteAccuracy Absolute accuracy. 079 * @param functionValueAccuracy Function value accuracy. 080 */ 081 public LaguerreSolver(double relativeAccuracy, 082 double absoluteAccuracy, 083 double functionValueAccuracy) { 084 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); 085 } 086 087 /** 088 * {@inheritDoc} 089 */ 090 @Override 091 public double doSolve() 092 throws TooManyEvaluationsException, 093 NumberIsTooLargeException, 094 NoBracketingException { 095 final double min = getMin(); 096 final double max = getMax(); 097 final double initial = getStartValue(); 098 final double functionValueAccuracy = getFunctionValueAccuracy(); 099 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, yMin, yInitial); 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, yInitial, yMax); 128 } 129 130 throw new NoBracketingException(min, max, yMin, yMax); 131 } 132 133 /** 134 * Find a real root in the given interval. 135 * 136 * Despite the bracketing condition, the root returned by 137 * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may 138 * not be a real zero inside {@code [min, max]}. 139 * For example, <code> p(x) = x<sup>3</sup> + 1, </code> 140 * with {@code min = -2}, {@code max = 2}, {@code initial = 0}. 141 * When it occurs, this code calls 142 * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)} 143 * in order to obtain all roots and picks up one real root. 144 * 145 * @param lo Lower bound of the search interval. 146 * @param hi Higher bound of the search interval. 147 * @param fLo Function value at the lower bound of the search interval. 148 * @param fHi Function value at the higher bound of the search interval. 149 * @return the point at which the function value is zero. 150 * @deprecated This method should not be part of the public API: It will 151 * be made private in version 4.0. 152 */ 153 @Deprecated 154 public double laguerre(double lo, double hi, 155 double fLo, double fHi) { 156 final Complex c[] = ComplexUtils.convertToComplex(getCoefficients()); 157 158 final Complex initial = new Complex(0.5 * (lo + hi), 0); 159 final Complex z = complexSolver.solve(c, initial); 160 if (complexSolver.isRoot(lo, hi, z)) { 161 return z.getReal(); 162 } else { 163 double r = Double.NaN; 164 // Solve all roots and select the one we are seeking. 165 Complex[] root = complexSolver.solveAll(c, initial); 166 for (int i = 0; i < root.length; i++) { 167 if (complexSolver.isRoot(lo, hi, root[i])) { 168 r = root[i].getReal(); 169 break; 170 } 171 } 172 return r; 173 } 174 } 175 176 /** 177 * Find all complex roots for the polynomial with the given 178 * coefficients, starting from the given initial value. 179 * <p> 180 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> 181 * 182 * @param coefficients Polynomial coefficients. 183 * @param initial Start value. 184 * @return the full set of complex roots of the polynomial 185 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException 186 * if the maximum number of evaluations is exceeded when solving for one of the roots 187 * @throws NullArgumentException if the {@code coefficients} is 188 * {@code null}. 189 * @throws NoDataException if the {@code coefficients} array is empty. 190 * @since 3.1 191 */ 192 public Complex[] solveAllComplex(double[] coefficients, 193 double initial) 194 throws NullArgumentException, 195 NoDataException, 196 TooManyEvaluationsException { 197 return solveAllComplex(coefficients, initial, Integer.MAX_VALUE); 198 } 199 200 /** 201 * Find all complex roots for the polynomial with the given 202 * coefficients, starting from the given initial value. 203 * <p> 204 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> 205 * 206 * @param coefficients polynomial coefficients 207 * @param initial start value 208 * @param maxEval maximum number of evaluations 209 * @return the full set of complex roots of the polynomial 210 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException 211 * if the maximum number of evaluations is exceeded when solving for one of the roots 212 * @throws NullArgumentException if the {@code coefficients} is 213 * {@code null} 214 * @throws NoDataException if the {@code coefficients} array is empty 215 * @since 3.5 216 */ 217 public Complex[] solveAllComplex(double[] coefficients, 218 double initial, int maxEval) 219 throws NullArgumentException, 220 NoDataException, 221 TooManyEvaluationsException { 222 setup(maxEval, 223 new PolynomialFunction(coefficients), 224 Double.NEGATIVE_INFINITY, 225 Double.POSITIVE_INFINITY, 226 initial); 227 return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients), 228 new Complex(initial, 0d)); 229 } 230 231 /** 232 * Find a complex root for the polynomial with the given coefficients, 233 * starting from the given initial value. 234 * <p> 235 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> 236 * 237 * @param coefficients Polynomial coefficients. 238 * @param initial Start value. 239 * @return a complex root of the polynomial 240 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException 241 * if the maximum number of evaluations is exceeded. 242 * @throws NullArgumentException if the {@code coefficients} is 243 * {@code null}. 244 * @throws NoDataException if the {@code coefficients} array is empty. 245 * @since 3.1 246 */ 247 public Complex solveComplex(double[] coefficients, 248 double initial) 249 throws NullArgumentException, 250 NoDataException, 251 TooManyEvaluationsException { 252 return solveComplex(coefficients, initial, Integer.MAX_VALUE); 253 } 254 255 /** 256 * Find a complex root for the polynomial with the given coefficients, 257 * starting from the given initial value. 258 * <p> 259 * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> 260 * 261 * @param coefficients polynomial coefficients 262 * @param initial start value 263 * @param maxEval maximum number of evaluations 264 * @return a complex root of the polynomial 265 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException 266 * if the maximum number of evaluations is exceeded 267 * @throws NullArgumentException if the {@code coefficients} is 268 * {@code null} 269 * @throws NoDataException if the {@code coefficients} array is empty 270 * @since 3.1 271 */ 272 public Complex solveComplex(double[] coefficients, 273 double initial, int maxEval) 274 throws NullArgumentException, 275 NoDataException, 276 TooManyEvaluationsException { 277 setup(maxEval, 278 new PolynomialFunction(coefficients), 279 Double.NEGATIVE_INFINITY, 280 Double.POSITIVE_INFINITY, 281 initial); 282 return complexSolver.solve(ComplexUtils.convertToComplex(coefficients), 283 new Complex(initial, 0d)); 284 } 285 286 /** 287 * Class for searching all (complex) roots. 288 */ 289 private class ComplexSolver { 290 /** 291 * Check whether the given complex root is actually a real zero 292 * in the given interval, within the solver tolerance level. 293 * 294 * @param min Lower bound for the interval. 295 * @param max Upper bound for the interval. 296 * @param z Complex root. 297 * @return {@code true} if z is a real zero. 298 */ 299 public boolean isRoot(double min, double max, Complex z) { 300 if (isSequence(min, z.getReal(), max)) { 301 double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy()); 302 return (FastMath.abs(z.getImaginary()) <= tolerance) || 303 (z.abs() <= getFunctionValueAccuracy()); 304 } 305 return false; 306 } 307 308 /** 309 * Find all complex roots for the polynomial with the given 310 * coefficients, 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[] solveAll(Complex coefficients[], Complex initial) 322 throws NullArgumentException, 323 NoDataException, 324 TooManyEvaluationsException { 325 if (coefficients == null) { 326 throw new NullArgumentException(); 327 } 328 final int n = coefficients.length - 1; 329 if (n == 0) { 330 throw new NoDataException(LocalizedFormats.POLYNOMIAL); 331 } 332 // Coefficients for deflated polynomial. 333 final Complex c[] = new Complex[n + 1]; 334 for (int i = 0; i <= n; i++) { 335 c[i] = coefficients[i]; 336 } 337 338 // Solve individual roots successively. 339 final Complex root[] = new Complex[n]; 340 for (int i = 0; i < n; i++) { 341 final Complex subarray[] = new Complex[n - i + 1]; 342 System.arraycopy(c, 0, subarray, 0, subarray.length); 343 root[i] = solve(subarray, initial); 344 // Polynomial deflation using synthetic division. 345 Complex newc = c[n - i]; 346 Complex oldc = null; 347 for (int j = n - i - 1; j >= 0; j--) { 348 oldc = c[j]; 349 c[j] = newc; 350 newc = oldc.add(newc.multiply(root[i])); 351 } 352 } 353 354 return root; 355 } 356 357 /** 358 * Find a complex root for the polynomial with the given coefficients, 359 * starting from the given initial value. 360 * 361 * @param coefficients Polynomial coefficients. 362 * @param initial Start value. 363 * @return the point at which the function value is zero. 364 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException 365 * if the maximum number of evaluations is exceeded. 366 * @throws NullArgumentException if the {@code coefficients} is 367 * {@code null}. 368 * @throws NoDataException if the {@code coefficients} array is empty. 369 */ 370 public Complex solve(Complex coefficients[], Complex initial) 371 throws NullArgumentException, 372 NoDataException, 373 TooManyEvaluationsException { 374 if (coefficients == null) { 375 throw new NullArgumentException(); 376 } 377 378 final int n = coefficients.length - 1; 379 if (n == 0) { 380 throw new NoDataException(LocalizedFormats.POLYNOMIAL); 381 } 382 383 final double absoluteAccuracy = getAbsoluteAccuracy(); 384 final double relativeAccuracy = getRelativeAccuracy(); 385 final double functionValueAccuracy = getFunctionValueAccuracy(); 386 387 final Complex nC = new Complex(n, 0); 388 final Complex n1C = new Complex(n - 1, 0); 389 390 Complex z = initial; 391 Complex oldz = new Complex(Double.POSITIVE_INFINITY, 392 Double.POSITIVE_INFINITY); 393 while (true) { 394 // Compute pv (polynomial value), dv (derivative value), and 395 // d2v (second derivative value) simultaneously. 396 Complex pv = coefficients[n]; 397 Complex dv = Complex.ZERO; 398 Complex d2v = Complex.ZERO; 399 for (int j = n-1; j >= 0; j--) { 400 d2v = dv.add(z.multiply(d2v)); 401 dv = pv.add(z.multiply(dv)); 402 pv = coefficients[j].add(z.multiply(pv)); 403 } 404 d2v = d2v.multiply(new Complex(2.0, 0.0)); 405 406 // Check for convergence. 407 final double tolerance = FastMath.max(relativeAccuracy * z.abs(), 408 absoluteAccuracy); 409 if ((z.subtract(oldz)).abs() <= tolerance) { 410 return z; 411 } 412 if (pv.abs() <= functionValueAccuracy) { 413 return z; 414 } 415 416 // Now pv != 0, calculate the new approximation. 417 final Complex G = dv.divide(pv); 418 final Complex G2 = G.multiply(G); 419 final Complex H = G2.subtract(d2v.divide(pv)); 420 final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2)); 421 // Choose a denominator larger in magnitude. 422 final Complex deltaSqrt = delta.sqrt(); 423 final Complex dplus = G.add(deltaSqrt); 424 final Complex dminus = G.subtract(deltaSqrt); 425 final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus; 426 // Perturb z if denominator is zero, for instance, 427 // p(x) = x^3 + 1, z = 0. 428 if (denominator.equals(new Complex(0.0, 0.0))) { 429 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy)); 430 oldz = new Complex(Double.POSITIVE_INFINITY, 431 Double.POSITIVE_INFINITY); 432 } else { 433 oldz = z; 434 z = z.subtract(nC.divide(denominator)); 435 } 436 incrementEvaluationCount(); 437 } 438 } 439 } 440}