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 * Except for the initial [min, max], it does not require bracketing 035 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If a complex 036 * number arises in the computation, we simply use its modulus as a real 037 * approximation.</p> 038 * <p> 039 * Because the interval may not be bracketing, the bisection alternative is 040 * not applicable here. However in practice our treatment usually works 041 * well, especially near real zeroes where the imaginary part of the complex 042 * approximation is often negligible.</p> 043 * <p> 044 * The formulas here do not use divided differences directly.</p> 045 * 046 * @since 1.2 047 * @see MullerSolver 048 */ 049public class MullerSolver2 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 MullerSolver2() { 058 this(DEFAULT_ABSOLUTE_ACCURACY); 059 } 060 /** 061 * Construct a solver. 062 * 063 * @param absoluteAccuracy Absolute accuracy. 064 */ 065 public MullerSolver2(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 MullerSolver2(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 090 verifyInterval(min, max); 091 092 final double relativeAccuracy = getRelativeAccuracy(); 093 final double absoluteAccuracy = getAbsoluteAccuracy(); 094 final double functionValueAccuracy = getFunctionValueAccuracy(); 095 096 // x2 is the last root approximation 097 // x is the new approximation and new x2 for next round 098 // x0 < x1 < x2 does not hold here 099 100 double x0 = min; 101 double y0 = computeObjectiveValue(x0); 102 if (FastMath.abs(y0) < functionValueAccuracy) { 103 return x0; 104 } 105 double x1 = max; 106 double y1 = computeObjectiveValue(x1); 107 if (FastMath.abs(y1) < functionValueAccuracy) { 108 return x1; 109 } 110 111 if(y0 * y1 > 0) { 112 throw new NoBracketingException(x0, x1, y0, y1); 113 } 114 115 double x2 = 0.5 * (x0 + x1); 116 double y2 = computeObjectiveValue(x2); 117 118 double oldx = Double.POSITIVE_INFINITY; 119 while (true) { 120 // quadratic interpolation through x0, x1, x2 121 final double q = (x2 - x1) / (x1 - x0); 122 final double a = q * (y2 - (1 + q) * y1 + q * y0); 123 final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; 124 final double c = (1 + q) * y2; 125 final double delta = b * b - 4 * a * c; 126 double x; 127 final double denominator; 128 if (delta >= 0.0) { 129 // choose a denominator larger in magnitude 130 double dplus = b + FastMath.sqrt(delta); 131 double dminus = b - FastMath.sqrt(delta); 132 denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus; 133 } else { 134 // take the modulus of (B +/- FastMath.sqrt(delta)) 135 denominator = FastMath.sqrt(b * b - delta); 136 } 137 if (denominator != 0) { 138 x = x2 - 2.0 * c * (x2 - x1) / denominator; 139 // perturb x if it exactly coincides with x1 or x2 140 // the equality tests here are intentional 141 while (x == x1 || x == x2) { 142 x += absoluteAccuracy; 143 } 144 } else { 145 // extremely rare case, get a random number to skip it 146 x = min + FastMath.random() * (max - min); 147 oldx = Double.POSITIVE_INFINITY; 148 } 149 final double y = computeObjectiveValue(x); 150 151 // check for convergence 152 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); 153 if (FastMath.abs(x - oldx) <= tolerance || 154 FastMath.abs(y) <= functionValueAccuracy) { 155 return x; 156 } 157 158 // prepare the next iteration 159 x0 = x1; 160 y0 = y1; 161 x1 = x2; 162 y1 = y2; 163 x2 = x; 164 y2 = y; 165 oldx = x; 166 } 167 } 168}