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.MaxIterationsExceededException; 22 import org.apache.commons.math.analysis.UnivariateRealFunction; 23 import org.apache.commons.math.util.FastMath; 24 import org.apache.commons.math.util.MathUtils; 25 26 /** 27 * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html"> 28 * Muller's Method</a> for root finding of real univariate functions. For 29 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, 30 * chapter 3. 31 * <p> 32 * Muller's method applies to both real and complex functions, but here we 33 * restrict ourselves to real functions. Methods solve() and solve2() find 34 * real zeros, using different ways to bypass complex arithmetics.</p> 35 * 36 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $ 37 * @since 1.2 38 */ 39 public class MullerSolver extends UnivariateRealSolverImpl { 40 41 /** 42 * Construct a solver for the given function. 43 * 44 * @param f function to solve 45 * @deprecated as of 2.0 the function to solve is passed as an argument 46 * to the {@link #solve(UnivariateRealFunction, double, double)} or 47 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} 48 * method. 49 */ 50 @Deprecated 51 public MullerSolver(UnivariateRealFunction f) { 52 super(f, 100, 1E-6); 53 } 54 55 /** 56 * Construct a solver. 57 * @deprecated in 2.2 (to be removed in 3.0). 58 */ 59 @Deprecated 60 public MullerSolver() { 61 super(100, 1E-6); 62 } 63 64 /** {@inheritDoc} */ 65 @Deprecated 66 public double solve(final double min, final double max) 67 throws ConvergenceException, FunctionEvaluationException { 68 return solve(f, min, max); 69 } 70 71 /** {@inheritDoc} */ 72 @Deprecated 73 public double solve(final double min, final double max, final double initial) 74 throws ConvergenceException, FunctionEvaluationException { 75 return solve(f, min, max, initial); 76 } 77 78 /** 79 * Find a real root in the given interval with initial value. 80 * <p> 81 * Requires bracketing condition.</p> 82 * 83 * @param f the function to solve 84 * @param min the lower bound for the interval 85 * @param max the upper bound for the interval 86 * @param initial the start value to use 87 * @param maxEval Maximum number of evaluations. 88 * @return the point at which the function value is zero 89 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 90 * or the solver detects convergence problems otherwise 91 * @throws FunctionEvaluationException if an error occurs evaluating the function 92 * @throws IllegalArgumentException if any parameters are invalid 93 */ 94 @Override 95 public double solve(int maxEval, final UnivariateRealFunction f, 96 final double min, final double max, final double initial) 97 throws MaxIterationsExceededException, FunctionEvaluationException { 98 setMaximalIterationCount(maxEval); 99 return solve(f, min, max, initial); 100 } 101 102 /** 103 * Find a real root in the given interval with initial value. 104 * <p> 105 * Requires bracketing condition.</p> 106 * 107 * @param f the function to solve 108 * @param min the lower bound for the interval 109 * @param max the upper bound for the interval 110 * @param initial the start value to use 111 * @return the point at which the function value is zero 112 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 113 * or the solver detects convergence problems otherwise 114 * @throws FunctionEvaluationException if an error occurs evaluating the function 115 * @throws IllegalArgumentException if any parameters are invalid 116 * @deprecated in 2.2 (to be removed in 3.0). 117 */ 118 @Deprecated 119 public double solve(final UnivariateRealFunction f, 120 final double min, final double max, final double initial) 121 throws MaxIterationsExceededException, FunctionEvaluationException { 122 123 // check for zeros before verifying bracketing 124 if (f.value(min) == 0.0) { return min; } 125 if (f.value(max) == 0.0) { return max; } 126 if (f.value(initial) == 0.0) { return initial; } 127 128 verifyBracketing(min, max, f); 129 verifySequence(min, initial, max); 130 if (isBracketing(min, initial, f)) { 131 return solve(f, min, initial); 132 } else { 133 return solve(f, initial, max); 134 } 135 } 136 137 /** 138 * Find a real root in the given interval. 139 * <p> 140 * Original Muller's method would have function evaluation at complex point. 141 * Since our f(x) is real, we have to find ways to avoid that. Bracketing 142 * condition is one way to go: by requiring bracketing in every iteration, 143 * the newly computed approximation is guaranteed to be real.</p> 144 * <p> 145 * Normally Muller's method converges quadratically in the vicinity of a 146 * zero, however it may be very slow in regions far away from zeros. For 147 * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use 148 * bisection as a safety backup if it performs very poorly.</p> 149 * <p> 150 * The formulas here use divided differences directly.</p> 151 * 152 * @param f the function to solve 153 * @param min the lower bound for the interval 154 * @param max the upper bound for the interval 155 * @param maxEval Maximum number of evaluations. 156 * @return the point at which the function value is zero 157 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 158 * or the solver detects convergence problems otherwise 159 * @throws FunctionEvaluationException if an error occurs evaluating the function 160 * @throws IllegalArgumentException if any parameters are invalid 161 */ 162 @Override 163 public double solve(int maxEval, final UnivariateRealFunction f, 164 final double min, final double max) 165 throws MaxIterationsExceededException, FunctionEvaluationException { 166 setMaximalIterationCount(maxEval); 167 return solve(f, min, max); 168 } 169 170 /** 171 * Find a real root in the given interval. 172 * <p> 173 * Original Muller's method would have function evaluation at complex point. 174 * Since our f(x) is real, we have to find ways to avoid that. Bracketing 175 * condition is one way to go: by requiring bracketing in every iteration, 176 * the newly computed approximation is guaranteed to be real.</p> 177 * <p> 178 * Normally Muller's method converges quadratically in the vicinity of a 179 * zero, however it may be very slow in regions far away from zeros. For 180 * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use 181 * bisection as a safety backup if it performs very poorly.</p> 182 * <p> 183 * The formulas here use divided differences directly.</p> 184 * 185 * @param f the function to solve 186 * @param min the lower bound for the interval 187 * @param max the upper bound for the interval 188 * @return the point at which the function value is zero 189 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 190 * or the solver detects convergence problems otherwise 191 * @throws FunctionEvaluationException if an error occurs evaluating the function 192 * @throws IllegalArgumentException if any parameters are invalid 193 * @deprecated in 2.2 (to be removed in 3.0). 194 */ 195 @Deprecated 196 public double solve(final UnivariateRealFunction f, 197 final double min, final double max) 198 throws MaxIterationsExceededException, FunctionEvaluationException { 199 200 // [x0, x2] is the bracketing interval in each iteration 201 // x1 is the last approximation and an interpolation point in (x0, x2) 202 // x is the new root approximation and new x1 for next round 203 // d01, d12, d012 are divided differences 204 205 double x0 = min; 206 double y0 = f.value(x0); 207 double x2 = max; 208 double y2 = f.value(x2); 209 double x1 = 0.5 * (x0 + x2); 210 double y1 = f.value(x1); 211 212 // check for zeros before verifying bracketing 213 if (y0 == 0.0) { 214 return min; 215 } 216 if (y2 == 0.0) { 217 return max; 218 } 219 verifyBracketing(min, max, f); 220 221 double oldx = Double.POSITIVE_INFINITY; 222 for (int i = 1; i <= maximalIterationCount; ++i) { 223 // Muller's method employs quadratic interpolation through 224 // x0, x1, x2 and x is the zero of the interpolating parabola. 225 // Due to bracketing condition, this parabola must have two 226 // real roots and we choose one in [x0, x2] to be x. 227 final double d01 = (y1 - y0) / (x1 - x0); 228 final double d12 = (y2 - y1) / (x2 - x1); 229 final double d012 = (d12 - d01) / (x2 - x0); 230 final double c1 = d01 + (x1 - x0) * d012; 231 final double delta = c1 * c1 - 4 * y1 * d012; 232 final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta)); 233 final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta)); 234 // xplus and xminus are two roots of parabola and at least 235 // one of them should lie in (x0, x2) 236 final double x = isSequence(x0, xplus, x2) ? xplus : xminus; 237 final double y = f.value(x); 238 239 // check for convergence 240 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); 241 if (FastMath.abs(x - oldx) <= tolerance) { 242 setResult(x, i); 243 return result; 244 } 245 if (FastMath.abs(y) <= functionValueAccuracy) { 246 setResult(x, i); 247 return result; 248 } 249 250 // Bisect if convergence is too slow. Bisection would waste 251 // our calculation of x, hopefully it won't happen often. 252 // the real number equality test x == x1 is intentional and 253 // completes the proximity tests above it 254 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || 255 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || 256 (x == x1); 257 // prepare the new bracketing interval for next iteration 258 if (!bisect) { 259 x0 = x < x1 ? x0 : x1; 260 y0 = x < x1 ? y0 : y1; 261 x2 = x > x1 ? x2 : x1; 262 y2 = x > x1 ? y2 : y1; 263 x1 = x; y1 = y; 264 oldx = x; 265 } else { 266 double xm = 0.5 * (x0 + x2); 267 double ym = f.value(xm); 268 if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) { 269 x2 = xm; y2 = ym; 270 } else { 271 x0 = xm; y0 = ym; 272 } 273 x1 = 0.5 * (x0 + x2); 274 y1 = f.value(x1); 275 oldx = Double.POSITIVE_INFINITY; 276 } 277 } 278 throw new MaxIterationsExceededException(maximalIterationCount); 279 } 280 281 /** 282 * Find a real root in the given interval. 283 * <p> 284 * solve2() differs from solve() in the way it avoids complex operations. 285 * Except for the initial [min, max], solve2() does not require bracketing 286 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex 287 * number arises in the computation, we simply use its modulus as real 288 * approximation.</p> 289 * <p> 290 * Because the interval may not be bracketing, bisection alternative is 291 * not applicable here. However in practice our treatment usually works 292 * well, especially near real zeros where the imaginary part of complex 293 * approximation is often negligible.</p> 294 * <p> 295 * The formulas here do not use divided differences directly.</p> 296 * 297 * @param min the lower bound for the interval 298 * @param max the upper bound for the interval 299 * @return the point at which the function value is zero 300 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 301 * or the solver detects convergence problems otherwise 302 * @throws FunctionEvaluationException if an error occurs evaluating the function 303 * @throws IllegalArgumentException if any parameters are invalid 304 * @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)} 305 * since 2.0 306 */ 307 @Deprecated 308 public double solve2(final double min, final double max) 309 throws MaxIterationsExceededException, FunctionEvaluationException { 310 return solve2(f, min, max); 311 } 312 313 /** 314 * Find a real root in the given interval. 315 * <p> 316 * solve2() differs from solve() in the way it avoids complex operations. 317 * Except for the initial [min, max], solve2() does not require bracketing 318 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex 319 * number arises in the computation, we simply use its modulus as real 320 * approximation.</p> 321 * <p> 322 * Because the interval may not be bracketing, bisection alternative is 323 * not applicable here. However in practice our treatment usually works 324 * well, especially near real zeros where the imaginary part of complex 325 * approximation is often negligible.</p> 326 * <p> 327 * The formulas here do not use divided differences directly.</p> 328 * 329 * @param f the function to solve 330 * @param min the lower bound for the interval 331 * @param max the upper bound for the interval 332 * @return the point at which the function value is zero 333 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded 334 * or the solver detects convergence problems otherwise 335 * @throws FunctionEvaluationException if an error occurs evaluating the function 336 * @throws IllegalArgumentException if any parameters are invalid 337 * @deprecated in 2.2 (to be removed in 3.0). 338 */ 339 @Deprecated 340 public double solve2(final UnivariateRealFunction f, 341 final double min, final double max) 342 throws MaxIterationsExceededException, FunctionEvaluationException { 343 344 // x2 is the last root approximation 345 // x is the new approximation and new x2 for next round 346 // x0 < x1 < x2 does not hold here 347 348 double x0 = min; 349 double y0 = f.value(x0); 350 double x1 = max; 351 double y1 = f.value(x1); 352 double x2 = 0.5 * (x0 + x1); 353 double y2 = f.value(x2); 354 355 // check for zeros before verifying bracketing 356 if (y0 == 0.0) { return min; } 357 if (y1 == 0.0) { return max; } 358 verifyBracketing(min, max, f); 359 360 double oldx = Double.POSITIVE_INFINITY; 361 for (int i = 1; i <= maximalIterationCount; ++i) { 362 // quadratic interpolation through x0, x1, x2 363 final double q = (x2 - x1) / (x1 - x0); 364 final double a = q * (y2 - (1 + q) * y1 + q * y0); 365 final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; 366 final double c = (1 + q) * y2; 367 final double delta = b * b - 4 * a * c; 368 double x; 369 final double denominator; 370 if (delta >= 0.0) { 371 // choose a denominator larger in magnitude 372 double dplus = b + FastMath.sqrt(delta); 373 double dminus = b - FastMath.sqrt(delta); 374 denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus; 375 } else { 376 // take the modulus of (B +/- FastMath.sqrt(delta)) 377 denominator = FastMath.sqrt(b * b - delta); 378 } 379 if (denominator != 0) { 380 x = x2 - 2.0 * c * (x2 - x1) / denominator; 381 // perturb x if it exactly coincides with x1 or x2 382 // the equality tests here are intentional 383 while (x == x1 || x == x2) { 384 x += absoluteAccuracy; 385 } 386 } else { 387 // extremely rare case, get a random number to skip it 388 x = min + FastMath.random() * (max - min); 389 oldx = Double.POSITIVE_INFINITY; 390 } 391 final double y = f.value(x); 392 393 // check for convergence 394 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); 395 if (FastMath.abs(x - oldx) <= tolerance) { 396 setResult(x, i); 397 return result; 398 } 399 if (FastMath.abs(y) <= functionValueAccuracy) { 400 setResult(x, i); 401 return result; 402 } 403 404 // prepare the next iteration 405 x0 = x1; 406 y0 = y1; 407 x1 = x2; 408 y1 = y2; 409 x2 = x; 410 y2 = y; 411 oldx = x; 412 } 413 throw new MaxIterationsExceededException(maximalIterationCount); 414 } 415 } 416