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