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