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 027 /** 028 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html"> 029 * Laguerre's Method</a> for root finding of real coefficient polynomials. 030 * For reference, see <b>A First Course in Numerical Analysis</b>, 031 * ISBN 048641454X, chapter 8. 032 * <p> 033 * Laguerre's method is global in the sense that it can start with any initial 034 * approximation and be able to solve all roots from that point.</p> 035 * 036 * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $ 037 * @since 1.2 038 */ 039 public class LaguerreSolver extends UnivariateRealSolverImpl { 040 041 /** Message for non-polynomial function. */ 042 private static final String NON_POLYNOMIAL_FUNCTION_MESSAGE = 043 "function is not polynomial"; 044 045 /** Message for non-positive degree. */ 046 private static final String NON_POSITIVE_DEGREE_MESSAGE = 047 "polynomial degree must be positive: degree={0}"; 048 049 /** polynomial function to solve. 050 * @deprecated as of 2.0 the function is not stored anymore in the instance 051 */ 052 @Deprecated 053 private final PolynomialFunction p; 054 055 /** 056 * Construct a solver for the given function. 057 * 058 * @param f function to solve 059 * @throws IllegalArgumentException if function is not polynomial 060 * @deprecated as of 2.0 the function to solve is passed as an argument 061 * to the {@link #solve(UnivariateRealFunction, double, double)} or 062 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 063 * method. 064 */ 065 @Deprecated 066 public LaguerreSolver(UnivariateRealFunction f) throws 067 IllegalArgumentException { 068 super(f, 100, 1E-6); 069 if (f instanceof PolynomialFunction) { 070 p = (PolynomialFunction) f; 071 } else { 072 throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE); 073 } 074 } 075 076 /** 077 * Construct a solver. 078 */ 079 public LaguerreSolver() { 080 super(100, 1E-6); 081 p = null; 082 } 083 084 /** 085 * Returns a copy of the polynomial function. 086 * 087 * @return a fresh copy of the polynomial function 088 * @deprecated as of 2.0 the function is not stored anymore within the instance. 089 */ 090 @Deprecated 091 public PolynomialFunction getPolynomialFunction() { 092 return new PolynomialFunction(p.getCoefficients()); 093 } 094 095 /** {@inheritDoc} */ 096 @Deprecated 097 public double solve(final double min, final double max) 098 throws ConvergenceException, FunctionEvaluationException { 099 return solve(p, min, max); 100 } 101 102 /** {@inheritDoc} */ 103 @Deprecated 104 public double solve(final double min, final double max, final double initial) 105 throws ConvergenceException, FunctionEvaluationException { 106 return solve(p, min, max, initial); 107 } 108 109 /** 110 * Find a real root in the given interval with initial value. 111 * <p> 112 * Requires bracketing condition.</p> 113 * 114 * @param f function to solve (must be polynomial) 115 * @param min the lower bound for the interval 116 * @param max the upper bound for the interval 117 * @param initial the start value to use 118 * @return the point at which the function value is zero 119 * @throws ConvergenceException if the maximum iteration count is exceeded 120 * or the solver detects convergence problems otherwise 121 * @throws FunctionEvaluationException if an error occurs evaluating the 122 * function 123 * @throws IllegalArgumentException if any parameters are invalid 124 */ 125 public double solve(final UnivariateRealFunction f, 126 final double min, final double max, final double initial) 127 throws ConvergenceException, FunctionEvaluationException { 128 129 // check for zeros before verifying bracketing 130 if (f.value(min) == 0.0) { 131 return min; 132 } 133 if (f.value(max) == 0.0) { 134 return max; 135 } 136 if (f.value(initial) == 0.0) { 137 return initial; 138 } 139 140 verifyBracketing(min, max, f); 141 verifySequence(min, initial, max); 142 if (isBracketing(min, initial, f)) { 143 return solve(f, min, initial); 144 } else { 145 return solve(f, initial, max); 146 } 147 148 } 149 150 /** 151 * Find a real root in the given interval. 152 * <p> 153 * Despite the bracketing condition, the root returned by solve(Complex[], 154 * Complex) may not be a real zero inside [min, max]. For example, 155 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try 156 * another initial value, or, as we did here, call solveAll() to obtain 157 * all roots and pick up the one that we're looking for.</p> 158 * 159 * @param f the function to solve 160 * @param min the lower bound for the interval 161 * @param max the upper bound for the interval 162 * @return the point at which the function value is zero 163 * @throws ConvergenceException if the maximum iteration count is exceeded 164 * or the solver detects convergence problems otherwise 165 * @throws FunctionEvaluationException if an error occurs evaluating the 166 * function 167 * @throws IllegalArgumentException if any parameters are invalid 168 */ 169 public double solve(final UnivariateRealFunction f, 170 final double min, final double max) 171 throws ConvergenceException, FunctionEvaluationException { 172 173 // check function type 174 if (!(f instanceof PolynomialFunction)) { 175 throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE); 176 } 177 178 // check for zeros before verifying bracketing 179 if (f.value(min) == 0.0) { return min; } 180 if (f.value(max) == 0.0) { return max; } 181 verifyBracketing(min, max, f); 182 183 double coefficients[] = ((PolynomialFunction) f).getCoefficients(); 184 Complex c[] = new Complex[coefficients.length]; 185 for (int i = 0; i < coefficients.length; i++) { 186 c[i] = new Complex(coefficients[i], 0.0); 187 } 188 Complex initial = new Complex(0.5 * (min + max), 0.0); 189 Complex z = solve(c, initial); 190 if (isRootOK(min, max, z)) { 191 setResult(z.getReal(), iterationCount); 192 return result; 193 } 194 195 // solve all roots and select the one we're seeking 196 Complex[] root = solveAll(c, initial); 197 for (int i = 0; i < root.length; i++) { 198 if (isRootOK(min, max, root[i])) { 199 setResult(root[i].getReal(), iterationCount); 200 return result; 201 } 202 } 203 204 // should never happen 205 throw new ConvergenceException(); 206 } 207 208 /** 209 * Returns true iff the given complex root is actually a real zero 210 * in the given interval, within the solver tolerance level. 211 * 212 * @param min the lower bound for the interval 213 * @param max the upper bound for the interval 214 * @param z the complex root 215 * @return true iff z is the sought-after real zero 216 */ 217 protected boolean isRootOK(double min, double max, Complex z) { 218 double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy); 219 return (isSequence(min, z.getReal(), max)) && 220 (Math.abs(z.getImaginary()) <= tolerance || 221 z.abs() <= functionValueAccuracy); 222 } 223 224 /** 225 * Find all complex roots for the polynomial with the given coefficients, 226 * starting from the given initial value. 227 * 228 * @param coefficients the polynomial coefficients array 229 * @param initial the start value to use 230 * @return the point at which the function value is zero 231 * @throws ConvergenceException if the maximum iteration count is exceeded 232 * or the solver detects convergence problems otherwise 233 * @throws FunctionEvaluationException if an error occurs evaluating the 234 * function 235 * @throws IllegalArgumentException if any parameters are invalid 236 */ 237 public Complex[] solveAll(double coefficients[], double initial) throws 238 ConvergenceException, FunctionEvaluationException { 239 240 Complex c[] = new Complex[coefficients.length]; 241 Complex z = new Complex(initial, 0.0); 242 for (int i = 0; i < c.length; i++) { 243 c[i] = new Complex(coefficients[i], 0.0); 244 } 245 return solveAll(c, z); 246 } 247 248 /** 249 * Find all complex roots for the polynomial with the given coefficients, 250 * starting from the given initial value. 251 * 252 * @param coefficients the polynomial coefficients array 253 * @param initial the start value to use 254 * @return the point at which the function value is zero 255 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 256 * or the solver detects convergence problems otherwise 257 * @throws FunctionEvaluationException if an error occurs evaluating the 258 * function 259 * @throws IllegalArgumentException if any parameters are invalid 260 */ 261 public Complex[] solveAll(Complex coefficients[], Complex initial) throws 262 MaxIterationsExceededException, FunctionEvaluationException { 263 264 int n = coefficients.length - 1; 265 int iterationCount = 0; 266 if (n < 1) { 267 throw MathRuntimeException.createIllegalArgumentException( 268 NON_POSITIVE_DEGREE_MESSAGE, n); 269 } 270 Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial 271 for (int i = 0; i <= n; i++) { 272 c[i] = coefficients[i]; 273 } 274 275 // solve individual root successively 276 Complex root[] = new Complex[n]; 277 for (int i = 0; i < n; i++) { 278 Complex subarray[] = new Complex[n-i+1]; 279 System.arraycopy(c, 0, subarray, 0, subarray.length); 280 root[i] = solve(subarray, initial); 281 // polynomial deflation using synthetic division 282 Complex newc = c[n-i]; 283 Complex oldc = null; 284 for (int j = n-i-1; j >= 0; j--) { 285 oldc = c[j]; 286 c[j] = newc; 287 newc = oldc.add(newc.multiply(root[i])); 288 } 289 iterationCount += this.iterationCount; 290 } 291 292 resultComputed = true; 293 this.iterationCount = iterationCount; 294 return root; 295 } 296 297 /** 298 * Find a complex root 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 307 * function 308 * @throws IllegalArgumentException if any parameters are invalid 309 */ 310 public Complex solve(Complex coefficients[], Complex initial) throws 311 MaxIterationsExceededException, FunctionEvaluationException { 312 313 int n = coefficients.length - 1; 314 if (n < 1) { 315 throw MathRuntimeException.createIllegalArgumentException( 316 NON_POSITIVE_DEGREE_MESSAGE, n); 317 } 318 Complex N = new Complex(n, 0.0); 319 Complex N1 = new Complex(n - 1, 0.0); 320 321 int i = 1; 322 Complex pv = null; 323 Complex dv = null; 324 Complex d2v = null; 325 Complex G = null; 326 Complex G2 = null; 327 Complex H = null; 328 Complex delta = null; 329 Complex denominator = null; 330 Complex z = initial; 331 Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY); 332 while (i <= maximalIterationCount) { 333 // Compute pv (polynomial value), dv (derivative value), and 334 // d2v (second derivative value) simultaneously. 335 pv = coefficients[n]; 336 dv = Complex.ZERO; 337 d2v = Complex.ZERO; 338 for (int j = n-1; j >= 0; j--) { 339 d2v = dv.add(z.multiply(d2v)); 340 dv = pv.add(z.multiply(dv)); 341 pv = coefficients[j].add(z.multiply(pv)); 342 } 343 d2v = d2v.multiply(new Complex(2.0, 0.0)); 344 345 // check for convergence 346 double tolerance = Math.max(relativeAccuracy * z.abs(), 347 absoluteAccuracy); 348 if ((z.subtract(oldz)).abs() <= tolerance) { 349 resultComputed = true; 350 iterationCount = i; 351 return z; 352 } 353 if (pv.abs() <= functionValueAccuracy) { 354 resultComputed = true; 355 iterationCount = i; 356 return z; 357 } 358 359 // now pv != 0, calculate the new approximation 360 G = dv.divide(pv); 361 G2 = G.multiply(G); 362 H = G2.subtract(d2v.divide(pv)); 363 delta = N1.multiply((N.multiply(H)).subtract(G2)); 364 // choose a denominator larger in magnitude 365 Complex deltaSqrt = delta.sqrt(); 366 Complex dplus = G.add(deltaSqrt); 367 Complex dminus = G.subtract(deltaSqrt); 368 denominator = dplus.abs() > dminus.abs() ? dplus : dminus; 369 // Perturb z if denominator is zero, for instance, 370 // p(x) = x^3 + 1, z = 0. 371 if (denominator.equals(new Complex(0.0, 0.0))) { 372 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy)); 373 oldz = new Complex(Double.POSITIVE_INFINITY, 374 Double.POSITIVE_INFINITY); 375 } else { 376 oldz = z; 377 z = z.subtract(N.divide(denominator)); 378 } 379 i++; 380 } 381 throw new MaxIterationsExceededException(maximalIterationCount); 382 } 383 }