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