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.math3.analysis.solvers; 018 019 import org.apache.commons.math3.util.FastMath; 020 import org.apache.commons.math3.exception.NumberIsTooLargeException; 021 import org.apache.commons.math3.exception.NoBracketingException; 022 import org.apache.commons.math3.exception.TooManyEvaluationsException; 023 024 /** 025 * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html"> 026 * Muller's Method</a> for root finding of real univariate functions. For 027 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, 028 * chapter 3. 029 * <p> 030 * Muller's method applies to both real and complex functions, but here we 031 * restrict ourselves to real functions. 032 * This class differs from {@link MullerSolver} in the way it avoids complex 033 * operations.</p> 034 * Muller's original method would have function evaluation at complex point. 035 * Since our f(x) is real, we have to find ways to avoid that. Bracketing 036 * condition is one way to go: by requiring bracketing in every iteration, 037 * the newly computed approximation is guaranteed to be real.</p> 038 * <p> 039 * Normally Muller's method converges quadratically in the vicinity of a 040 * zero, however it may be very slow in regions far away from zeros. For 041 * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use 042 * bisection as a safety backup if it performs very poorly.</p> 043 * <p> 044 * The formulas here use divided differences directly.</p> 045 * 046 * @version $Id: MullerSolver.java 1391927 2012-09-30 00:03:30Z erans $ 047 * @since 1.2 048 * @see MullerSolver2 049 */ 050 public class MullerSolver extends AbstractUnivariateSolver { 051 052 /** Default absolute accuracy. */ 053 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; 054 055 /** 056 * Construct a solver with default accuracy (1e-6). 057 */ 058 public MullerSolver() { 059 this(DEFAULT_ABSOLUTE_ACCURACY); 060 } 061 /** 062 * Construct a solver. 063 * 064 * @param absoluteAccuracy Absolute accuracy. 065 */ 066 public MullerSolver(double absoluteAccuracy) { 067 super(absoluteAccuracy); 068 } 069 /** 070 * Construct a solver. 071 * 072 * @param relativeAccuracy Relative accuracy. 073 * @param absoluteAccuracy Absolute accuracy. 074 */ 075 public MullerSolver(double relativeAccuracy, 076 double absoluteAccuracy) { 077 super(relativeAccuracy, absoluteAccuracy); 078 } 079 080 /** 081 * {@inheritDoc} 082 */ 083 @Override 084 protected double doSolve() 085 throws TooManyEvaluationsException, 086 NumberIsTooLargeException, 087 NoBracketingException { 088 final double min = getMin(); 089 final double max = getMax(); 090 final double initial = getStartValue(); 091 092 final double functionValueAccuracy = getFunctionValueAccuracy(); 093 094 verifySequence(min, initial, max); 095 096 // check for zeros before verifying bracketing 097 final double fMin = computeObjectiveValue(min); 098 if (FastMath.abs(fMin) < functionValueAccuracy) { 099 return min; 100 } 101 final double fMax = computeObjectiveValue(max); 102 if (FastMath.abs(fMax) < functionValueAccuracy) { 103 return max; 104 } 105 final double fInitial = computeObjectiveValue(initial); 106 if (FastMath.abs(fInitial) < functionValueAccuracy) { 107 return initial; 108 } 109 110 verifyBracketing(min, max); 111 112 if (isBracketing(min, initial)) { 113 return solve(min, initial, fMin, fInitial); 114 } else { 115 return solve(initial, max, fInitial, fMax); 116 } 117 } 118 119 /** 120 * Find a real root in the given interval. 121 * 122 * @param min Lower bound for the interval. 123 * @param max Upper bound for the interval. 124 * @param fMin function value at the lower bound. 125 * @param fMax function value at the upper bound. 126 * @return the point at which the function value is zero. 127 * @throws TooManyEvaluationsException if the allowed number of calls to 128 * the function to be solved has been exhausted. 129 */ 130 private double solve(double min, double max, 131 double fMin, double fMax) 132 throws TooManyEvaluationsException { 133 final double relativeAccuracy = getRelativeAccuracy(); 134 final double absoluteAccuracy = getAbsoluteAccuracy(); 135 final double functionValueAccuracy = getFunctionValueAccuracy(); 136 137 // [x0, x2] is the bracketing interval in each iteration 138 // x1 is the last approximation and an interpolation point in (x0, x2) 139 // x is the new root approximation and new x1 for next round 140 // d01, d12, d012 are divided differences 141 142 double x0 = min; 143 double y0 = fMin; 144 double x2 = max; 145 double y2 = fMax; 146 double x1 = 0.5 * (x0 + x2); 147 double y1 = computeObjectiveValue(x1); 148 149 double oldx = Double.POSITIVE_INFINITY; 150 while (true) { 151 // Muller's method employs quadratic interpolation through 152 // x0, x1, x2 and x is the zero of the interpolating parabola. 153 // Due to bracketing condition, this parabola must have two 154 // real roots and we choose one in [x0, x2] to be x. 155 final double d01 = (y1 - y0) / (x1 - x0); 156 final double d12 = (y2 - y1) / (x2 - x1); 157 final double d012 = (d12 - d01) / (x2 - x0); 158 final double c1 = d01 + (x1 - x0) * d012; 159 final double delta = c1 * c1 - 4 * y1 * d012; 160 final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta)); 161 final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta)); 162 // xplus and xminus are two roots of parabola and at least 163 // one of them should lie in (x0, x2) 164 final double x = isSequence(x0, xplus, x2) ? xplus : xminus; 165 final double y = computeObjectiveValue(x); 166 167 // check for convergence 168 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); 169 if (FastMath.abs(x - oldx) <= tolerance || 170 FastMath.abs(y) <= functionValueAccuracy) { 171 return x; 172 } 173 174 // Bisect if convergence is too slow. Bisection would waste 175 // our calculation of x, hopefully it won't happen often. 176 // the real number equality test x == x1 is intentional and 177 // completes the proximity tests above it 178 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || 179 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || 180 (x == x1); 181 // prepare the new bracketing interval for next iteration 182 if (!bisect) { 183 x0 = x < x1 ? x0 : x1; 184 y0 = x < x1 ? y0 : y1; 185 x2 = x > x1 ? x2 : x1; 186 y2 = x > x1 ? y2 : y1; 187 x1 = x; y1 = y; 188 oldx = x; 189 } else { 190 double xm = 0.5 * (x0 + x2); 191 double ym = computeObjectiveValue(xm); 192 if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) { 193 x2 = xm; y2 = ym; 194 } else { 195 x0 = xm; y0 = ym; 196 } 197 x1 = 0.5 * (x0 + x2); 198 y1 = computeObjectiveValue(x1); 199 oldx = Double.POSITIVE_INFINITY; 200 } 201 } 202 } 203 }