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 */ 017package org.apache.commons.math3.analysis.solvers; 018 019import org.apache.commons.math3.exception.NoBracketingException; 020import org.apache.commons.math3.exception.NumberIsTooLargeException; 021import org.apache.commons.math3.exception.TooManyEvaluationsException; 022import org.apache.commons.math3.util.FastMath; 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><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 * @since 1.2 047 * @see MullerSolver2 048 */ 049public class MullerSolver extends AbstractUnivariateSolver { 050 051 /** Default absolute accuracy. */ 052 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; 053 054 /** 055 * Construct a solver with default accuracy (1e-6). 056 */ 057 public MullerSolver() { 058 this(DEFAULT_ABSOLUTE_ACCURACY); 059 } 060 /** 061 * Construct a solver. 062 * 063 * @param absoluteAccuracy Absolute accuracy. 064 */ 065 public MullerSolver(double absoluteAccuracy) { 066 super(absoluteAccuracy); 067 } 068 /** 069 * Construct a solver. 070 * 071 * @param relativeAccuracy Relative accuracy. 072 * @param absoluteAccuracy Absolute accuracy. 073 */ 074 public MullerSolver(double relativeAccuracy, 075 double absoluteAccuracy) { 076 super(relativeAccuracy, absoluteAccuracy); 077 } 078 079 /** 080 * {@inheritDoc} 081 */ 082 @Override 083 protected double doSolve() 084 throws TooManyEvaluationsException, 085 NumberIsTooLargeException, 086 NoBracketingException { 087 final double min = getMin(); 088 final double max = getMax(); 089 final double initial = getStartValue(); 090 091 final double functionValueAccuracy = getFunctionValueAccuracy(); 092 093 verifySequence(min, initial, max); 094 095 // check for zeros before verifying bracketing 096 final double fMin = computeObjectiveValue(min); 097 if (FastMath.abs(fMin) < functionValueAccuracy) { 098 return min; 099 } 100 final double fMax = computeObjectiveValue(max); 101 if (FastMath.abs(fMax) < functionValueAccuracy) { 102 return max; 103 } 104 final double fInitial = computeObjectiveValue(initial); 105 if (FastMath.abs(fInitial) < functionValueAccuracy) { 106 return initial; 107 } 108 109 verifyBracketing(min, max); 110 111 if (isBracketing(min, initial)) { 112 return solve(min, initial, fMin, fInitial); 113 } else { 114 return solve(initial, max, fInitial, fMax); 115 } 116 } 117 118 /** 119 * Find a real root in the given interval. 120 * 121 * @param min Lower bound for the interval. 122 * @param max Upper bound for the interval. 123 * @param fMin function value at the lower bound. 124 * @param fMax function value at the upper bound. 125 * @return the point at which the function value is zero. 126 * @throws TooManyEvaluationsException if the allowed number of calls to 127 * the function to be solved has been exhausted. 128 */ 129 private double solve(double min, double max, 130 double fMin, double fMax) 131 throws TooManyEvaluationsException { 132 final double relativeAccuracy = getRelativeAccuracy(); 133 final double absoluteAccuracy = getAbsoluteAccuracy(); 134 final double functionValueAccuracy = getFunctionValueAccuracy(); 135 136 // [x0, x2] is the bracketing interval in each iteration 137 // x1 is the last approximation and an interpolation point in (x0, x2) 138 // x is the new root approximation and new x1 for next round 139 // d01, d12, d012 are divided differences 140 141 double x0 = min; 142 double y0 = fMin; 143 double x2 = max; 144 double y2 = fMax; 145 double x1 = 0.5 * (x0 + x2); 146 double y1 = computeObjectiveValue(x1); 147 148 double oldx = Double.POSITIVE_INFINITY; 149 while (true) { 150 // Muller's method employs quadratic interpolation through 151 // x0, x1, x2 and x is the zero of the interpolating parabola. 152 // Due to bracketing condition, this parabola must have two 153 // real roots and we choose one in [x0, x2] to be x. 154 final double d01 = (y1 - y0) / (x1 - x0); 155 final double d12 = (y2 - y1) / (x2 - x1); 156 final double d012 = (d12 - d01) / (x2 - x0); 157 final double c1 = d01 + (x1 - x0) * d012; 158 final double delta = c1 * c1 - 4 * y1 * d012; 159 final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta)); 160 final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta)); 161 // xplus and xminus are two roots of parabola and at least 162 // one of them should lie in (x0, x2) 163 final double x = isSequence(x0, xplus, x2) ? xplus : xminus; 164 final double y = computeObjectiveValue(x); 165 166 // check for convergence 167 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); 168 if (FastMath.abs(x - oldx) <= tolerance || 169 FastMath.abs(y) <= functionValueAccuracy) { 170 return x; 171 } 172 173 // Bisect if convergence is too slow. Bisection would waste 174 // our calculation of x, hopefully it won't happen often. 175 // the real number equality test x == x1 is intentional and 176 // completes the proximity tests above it 177 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || 178 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || 179 (x == x1); 180 // prepare the new bracketing interval for next iteration 181 if (!bisect) { 182 x0 = x < x1 ? x0 : x1; 183 y0 = x < x1 ? y0 : y1; 184 x2 = x > x1 ? x2 : x1; 185 y2 = x > x1 ? y2 : y1; 186 x1 = x; y1 = y; 187 oldx = x; 188 } else { 189 double xm = 0.5 * (x0 + x2); 190 double ym = computeObjectiveValue(xm); 191 if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) { 192 x2 = xm; y2 = ym; 193 } else { 194 x0 = xm; y0 = ym; 195 } 196 x1 = 0.5 * (x0 + x2); 197 y1 = computeObjectiveValue(x1); 198 oldx = Double.POSITIVE_INFINITY; 199 } 200 } 201 } 202}