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 * http://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 package org.apache.commons.math.analysis.solvers; 18 19 import org.apache.commons.math.ConvergenceException; 20 import org.apache.commons.math.FunctionEvaluationException; 21 import org.apache.commons.math.MathRuntimeException; 22 import org.apache.commons.math.MaxIterationsExceededException; 23 import org.apache.commons.math.analysis.UnivariateRealFunction; 24 import org.apache.commons.math.analysis.polynomials.PolynomialFunction; 25 import org.apache.commons.math.complex.Complex; 26 import org.apache.commons.math.exception.util.LocalizedFormats; 27 import org.apache.commons.math.util.FastMath; 28 29 /** 30 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html"> 31 * Laguerre's Method</a> for root finding of real coefficient polynomials. 32 * For reference, see <b>A First Course in Numerical Analysis</b>, 33 * ISBN 048641454X, chapter 8. 34 * <p> 35 * Laguerre's method is global in the sense that it can start with any initial 36 * approximation and be able to solve all roots from that point.</p> 37 * 38 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $ 39 * @since 1.2 40 */ 41 public class LaguerreSolver extends UnivariateRealSolverImpl { 42 43 /** polynomial function to solve. 44 * @deprecated as of 2.0 the function is not stored anymore in the instance 45 */ 46 @Deprecated 47 private final PolynomialFunction p; 48 49 /** 50 * Construct a solver for the given function. 51 * 52 * @param f function to solve 53 * @throws IllegalArgumentException if function is not polynomial 54 * @deprecated as of 2.0 the function to solve is passed as an argument 55 * to the {@link #solve(UnivariateRealFunction, double, double)} or 56 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 57 * method. 58 */ 59 @Deprecated 60 public LaguerreSolver(UnivariateRealFunction f) throws IllegalArgumentException { 61 super(f, 100, 1E-6); 62 if (f instanceof PolynomialFunction) { 63 p = (PolynomialFunction) f; 64 } else { 65 throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL); 66 } 67 } 68 69 /** 70 * Construct a solver. 71 * @deprecated in 2.2 (to be removed in 3.0) 72 */ 73 @Deprecated 74 public LaguerreSolver() { 75 super(100, 1E-6); 76 p = null; 77 } 78 79 /** 80 * Returns a copy of the polynomial function. 81 * 82 * @return a fresh copy of the polynomial function 83 * @deprecated as of 2.0 the function is not stored anymore within the instance. 84 */ 85 @Deprecated 86 public PolynomialFunction getPolynomialFunction() { 87 return new PolynomialFunction(p.getCoefficients()); 88 } 89 90 /** {@inheritDoc} */ 91 @Deprecated 92 public double solve(final double min, final double max) 93 throws ConvergenceException, FunctionEvaluationException { 94 return solve(p, min, max); 95 } 96 97 /** {@inheritDoc} */ 98 @Deprecated 99 public double solve(final double min, final double max, final double initial) 100 throws ConvergenceException, FunctionEvaluationException { 101 return solve(p, min, max, initial); 102 } 103 104 /** 105 * Find a real root in the given interval with initial value. 106 * <p> 107 * Requires bracketing condition.</p> 108 * 109 * @param f function to solve (must be polynomial) 110 * @param min the lower bound for the interval 111 * @param max the upper bound for the interval 112 * @param initial the start value to use 113 * @param maxEval Maximum number of evaluations. 114 * @return the point at which the function value is zero 115 * @throws ConvergenceException if the maximum iteration count is exceeded 116 * or the solver detects convergence problems otherwise 117 * @throws FunctionEvaluationException if an error occurs evaluating the function 118 * @throws IllegalArgumentException if any parameters are invalid 119 */ 120 @Override 121 public double solve(int maxEval, final UnivariateRealFunction f, 122 final double min, final double max, final double initial) 123 throws ConvergenceException, FunctionEvaluationException { 124 setMaximalIterationCount(maxEval); 125 return solve(f, min, max, initial); 126 } 127 128 /** 129 * Find a real root in the given interval with initial value. 130 * <p> 131 * Requires bracketing condition.</p> 132 * 133 * @param f function to solve (must be polynomial) 134 * @param min the lower bound for the interval 135 * @param max the upper bound for the interval 136 * @param initial the start value to use 137 * @return the point at which the function value is zero 138 * @throws ConvergenceException if the maximum iteration count is exceeded 139 * or the solver detects convergence problems otherwise 140 * @throws FunctionEvaluationException if an error occurs evaluating the function 141 * @throws IllegalArgumentException if any parameters are invalid 142 * @deprecated in 2.2 (to be removed in 3.0). 143 */ 144 @Deprecated 145 public double solve(final UnivariateRealFunction f, 146 final double min, final double max, final double initial) 147 throws ConvergenceException, FunctionEvaluationException { 148 149 // check for zeros before verifying bracketing 150 if (f.value(min) == 0.0) { 151 return min; 152 } 153 if (f.value(max) == 0.0) { 154 return max; 155 } 156 if (f.value(initial) == 0.0) { 157 return initial; 158 } 159 160 verifyBracketing(min, max, f); 161 verifySequence(min, initial, max); 162 if (isBracketing(min, initial, f)) { 163 return solve(f, min, initial); 164 } else { 165 return solve(f, initial, max); 166 } 167 168 } 169 170 /** 171 * Find a real root in the given interval. 172 * <p> 173 * Despite the bracketing condition, the root returned by solve(Complex[], 174 * Complex) may not be a real zero inside [min, max]. For example, 175 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try 176 * another initial value, or, as we did here, call solveAll() to obtain 177 * all roots and pick up the one that we're looking for.</p> 178 * 179 * @param f the function to solve 180 * @param min the lower bound for the interval 181 * @param max the upper bound for the interval 182 * @param maxEval Maximum number of evaluations. 183 * @return the point at which the function value is zero 184 * @throws ConvergenceException if the maximum iteration count is exceeded 185 * or the solver detects convergence problems otherwise 186 * @throws FunctionEvaluationException if an error occurs evaluating the function 187 * @throws IllegalArgumentException if any parameters are invalid 188 */ 189 @Override 190 public double solve(int maxEval, final UnivariateRealFunction f, 191 final double min, final double max) 192 throws ConvergenceException, FunctionEvaluationException { 193 setMaximalIterationCount(maxEval); 194 return solve(f, min, max); 195 } 196 197 /** 198 * Find a real root in the given interval. 199 * <p> 200 * Despite the bracketing condition, the root returned by solve(Complex[], 201 * Complex) may not be a real zero inside [min, max]. For example, 202 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try 203 * another initial value, or, as we did here, call solveAll() to obtain 204 * all roots and pick up the one that we're looking for.</p> 205 * 206 * @param f the function to solve 207 * @param min the lower bound for the interval 208 * @param max the upper bound for the interval 209 * @return the point at which the function value is zero 210 * @throws ConvergenceException if the maximum iteration count is exceeded 211 * or the solver detects convergence problems otherwise 212 * @throws FunctionEvaluationException if an error occurs evaluating the function 213 * @throws IllegalArgumentException if any parameters are invalid 214 * @deprecated in 2.2 (to be removed in 3.0). 215 */ 216 @Deprecated 217 public double solve(final UnivariateRealFunction f, 218 final double min, final double max) 219 throws ConvergenceException, FunctionEvaluationException { 220 221 // check function type 222 if (!(f instanceof PolynomialFunction)) { 223 throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL); 224 } 225 226 // check for zeros before verifying bracketing 227 if (f.value(min) == 0.0) { return min; } 228 if (f.value(max) == 0.0) { return max; } 229 verifyBracketing(min, max, f); 230 231 double coefficients[] = ((PolynomialFunction) f).getCoefficients(); 232 Complex c[] = new Complex[coefficients.length]; 233 for (int i = 0; i < coefficients.length; i++) { 234 c[i] = new Complex(coefficients[i], 0.0); 235 } 236 Complex initial = new Complex(0.5 * (min + max), 0.0); 237 Complex z = solve(c, initial); 238 if (isRootOK(min, max, z)) { 239 setResult(z.getReal(), iterationCount); 240 return result; 241 } 242 243 // solve all roots and select the one we're seeking 244 Complex[] root = solveAll(c, initial); 245 for (int i = 0; i < root.length; i++) { 246 if (isRootOK(min, max, root[i])) { 247 setResult(root[i].getReal(), iterationCount); 248 return result; 249 } 250 } 251 252 // should never happen 253 throw new ConvergenceException(); 254 } 255 256 /** 257 * Returns true iff the given complex root is actually a real zero 258 * in the given interval, within the solver tolerance level. 259 * 260 * @param min the lower bound for the interval 261 * @param max the upper bound for the interval 262 * @param z the complex root 263 * @return true iff z is the sought-after real zero 264 */ 265 protected boolean isRootOK(double min, double max, Complex z) { 266 double tolerance = FastMath.max(relativeAccuracy * z.abs(), absoluteAccuracy); 267 return (isSequence(min, z.getReal(), max)) && 268 (FastMath.abs(z.getImaginary()) <= tolerance || 269 z.abs() <= functionValueAccuracy); 270 } 271 272 /** 273 * Find all complex roots for the polynomial with the given coefficients, 274 * starting from the given initial value. 275 * 276 * @param coefficients the polynomial coefficients array 277 * @param initial the start value to use 278 * @return the point at which the function value is zero 279 * @throws ConvergenceException if the maximum iteration count is exceeded 280 * or the solver detects convergence problems otherwise 281 * @throws FunctionEvaluationException if an error occurs evaluating the function 282 * @throws IllegalArgumentException if any parameters are invalid 283 * @deprecated in 2.2. 284 */ 285 @Deprecated 286 public Complex[] solveAll(double coefficients[], double initial) throws 287 ConvergenceException, FunctionEvaluationException { 288 289 Complex c[] = new Complex[coefficients.length]; 290 Complex z = new Complex(initial, 0.0); 291 for (int i = 0; i < c.length; i++) { 292 c[i] = new Complex(coefficients[i], 0.0); 293 } 294 return solveAll(c, z); 295 } 296 297 /** 298 * Find all complex roots for the polynomial with the given coefficients, 299 * starting from the given initial value. 300 * 301 * @param coefficients the polynomial coefficients array 302 * @param initial the start value to use 303 * @return the point at which the function value is zero 304 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 305 * or the solver detects convergence problems otherwise 306 * @throws FunctionEvaluationException if an error occurs evaluating the function 307 * @throws IllegalArgumentException if any parameters are invalid 308 * @deprecated in 2.2. 309 */ 310 @Deprecated 311 public Complex[] solveAll(Complex coefficients[], Complex initial) throws 312 MaxIterationsExceededException, FunctionEvaluationException { 313 314 int n = coefficients.length - 1; 315 int iterationCount = 0; 316 if (n < 1) { 317 throw MathRuntimeException.createIllegalArgumentException( 318 LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n); 319 } 320 Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial 321 for (int i = 0; i <= n; i++) { 322 c[i] = coefficients[i]; 323 } 324 325 // solve individual root successively 326 Complex root[] = new Complex[n]; 327 for (int i = 0; i < n; i++) { 328 Complex subarray[] = new Complex[n-i+1]; 329 System.arraycopy(c, 0, subarray, 0, subarray.length); 330 root[i] = solve(subarray, initial); 331 // polynomial deflation using synthetic division 332 Complex newc = c[n-i]; 333 Complex oldc = null; 334 for (int j = n-i-1; j >= 0; j--) { 335 oldc = c[j]; 336 c[j] = newc; 337 newc = oldc.add(newc.multiply(root[i])); 338 } 339 iterationCount += this.iterationCount; 340 } 341 342 resultComputed = true; 343 this.iterationCount = iterationCount; 344 return root; 345 } 346 347 /** 348 * Find a complex root for the polynomial with the given coefficients, 349 * starting from the given initial value. 350 * 351 * @param coefficients the polynomial coefficients array 352 * @param initial the start value to use 353 * @return the point at which the function value is zero 354 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 355 * or the solver detects convergence problems otherwise 356 * @throws FunctionEvaluationException if an error occurs evaluating the function 357 * @throws IllegalArgumentException if any parameters are invalid 358 * @deprecated in 2.2. 359 */ 360 @Deprecated 361 public Complex solve(Complex coefficients[], Complex initial) throws 362 MaxIterationsExceededException, FunctionEvaluationException { 363 364 int n = coefficients.length - 1; 365 if (n < 1) { 366 throw MathRuntimeException.createIllegalArgumentException( 367 LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n); 368 } 369 Complex N = new Complex(n, 0.0); 370 Complex N1 = new Complex(n - 1, 0.0); 371 372 int i = 1; 373 Complex pv = null; 374 Complex dv = null; 375 Complex d2v = null; 376 Complex G = null; 377 Complex G2 = null; 378 Complex H = null; 379 Complex delta = null; 380 Complex denominator = null; 381 Complex z = initial; 382 Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY); 383 while (i <= maximalIterationCount) { 384 // Compute pv (polynomial value), dv (derivative value), and 385 // d2v (second derivative value) simultaneously. 386 pv = coefficients[n]; 387 dv = Complex.ZERO; 388 d2v = Complex.ZERO; 389 for (int j = n-1; j >= 0; j--) { 390 d2v = dv.add(z.multiply(d2v)); 391 dv = pv.add(z.multiply(dv)); 392 pv = coefficients[j].add(z.multiply(pv)); 393 } 394 d2v = d2v.multiply(new Complex(2.0, 0.0)); 395 396 // check for convergence 397 double tolerance = FastMath.max(relativeAccuracy * z.abs(), 398 absoluteAccuracy); 399 if ((z.subtract(oldz)).abs() <= tolerance) { 400 resultComputed = true; 401 iterationCount = i; 402 return z; 403 } 404 if (pv.abs() <= functionValueAccuracy) { 405 resultComputed = true; 406 iterationCount = i; 407 return z; 408 } 409 410 // now pv != 0, calculate the new approximation 411 G = dv.divide(pv); 412 G2 = G.multiply(G); 413 H = G2.subtract(d2v.divide(pv)); 414 delta = N1.multiply((N.multiply(H)).subtract(G2)); 415 // choose a denominator larger in magnitude 416 Complex deltaSqrt = delta.sqrt(); 417 Complex dplus = G.add(deltaSqrt); 418 Complex dminus = G.subtract(deltaSqrt); 419 denominator = dplus.abs() > dminus.abs() ? dplus : dminus; 420 // Perturb z if denominator is zero, for instance, 421 // p(x) = x^3 + 1, z = 0. 422 if (denominator.equals(new Complex(0.0, 0.0))) { 423 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy)); 424 oldz = new Complex(Double.POSITIVE_INFINITY, 425 Double.POSITIVE_INFINITY); 426 } else { 427 oldz = z; 428 z = z.subtract(N.divide(denominator)); 429 } 430 i++; 431 } 432 throw new MaxIterationsExceededException(maximalIterationCount); 433 } 434 } 435
