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.math.analysis.solvers; 018 019 import org.apache.commons.math.ConvergenceException; 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.MathRuntimeException; 022 import org.apache.commons.math.MaxIterationsExceededException; 023 import org.apache.commons.math.analysis.UnivariateRealFunction; 024 import org.apache.commons.math.analysis.polynomials.PolynomialFunction; 025 import org.apache.commons.math.complex.Complex; 026 import org.apache.commons.math.exception.util.LocalizedFormats; 027 import org.apache.commons.math.util.FastMath; 028 029 /** 030 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html"> 031 * Laguerre's Method</a> for root finding of real coefficient polynomials. 032 * For reference, see <b>A First Course in Numerical Analysis</b>, 033 * ISBN 048641454X, chapter 8. 034 * <p> 035 * Laguerre's method is global in the sense that it can start with any initial 036 * approximation and be able to solve all roots from that point.</p> 037 * 038 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $ 039 * @since 1.2 040 */ 041 public class LaguerreSolver extends UnivariateRealSolverImpl { 042 043 /** polynomial function to solve. 044 * @deprecated as of 2.0 the function is not stored anymore in the instance 045 */ 046 @Deprecated 047 private final PolynomialFunction p; 048 049 /** 050 * Construct a solver for the given function. 051 * 052 * @param f function to solve 053 * @throws IllegalArgumentException if function is not polynomial 054 * @deprecated as of 2.0 the function to solve is passed as an argument 055 * to the {@link #solve(UnivariateRealFunction, double, double)} or 056 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 057 * method. 058 */ 059 @Deprecated 060 public LaguerreSolver(UnivariateRealFunction f) throws IllegalArgumentException { 061 super(f, 100, 1E-6); 062 if (f instanceof PolynomialFunction) { 063 p = (PolynomialFunction) f; 064 } else { 065 throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL); 066 } 067 } 068 069 /** 070 * Construct a solver. 071 * @deprecated in 2.2 (to be removed in 3.0) 072 */ 073 @Deprecated 074 public LaguerreSolver() { 075 super(100, 1E-6); 076 p = null; 077 } 078 079 /** 080 * Returns a copy of the polynomial function. 081 * 082 * @return a fresh copy of the polynomial function 083 * @deprecated as of 2.0 the function is not stored anymore within the instance. 084 */ 085 @Deprecated 086 public PolynomialFunction getPolynomialFunction() { 087 return new PolynomialFunction(p.getCoefficients()); 088 } 089 090 /** {@inheritDoc} */ 091 @Deprecated 092 public double solve(final double min, final double max) 093 throws ConvergenceException, FunctionEvaluationException { 094 return solve(p, min, max); 095 } 096 097 /** {@inheritDoc} */ 098 @Deprecated 099 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 }