http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java deleted file mode 100644 index 4981592..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java +++ /dev/null @@ -1,412 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - - -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.exception.MathInternalError; -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.Precision; - -/** - * This class implements a modification of the <a - * href="http://mathworld.wolfram.com/BrentsMethod.html";> Brent algorithm</a>. - * <p> - * The changes with respect to the original Brent algorithm are: - * <ul> - * <li>the returned value is chosen in the current interval according - * to user specified {@link AllowedSolution},</li> - * <li>the maximal order for the invert polynomial root search is - * user-specified instead of being invert quadratic only</li> - * </ul> - * </p> - * The given interval must bracket the root. - * - */ -public class BracketingNthOrderBrentSolver - extends AbstractUnivariateSolver - implements BracketedUnivariateSolver<UnivariateFunction> { - - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** Default maximal order. */ - private static final int DEFAULT_MAXIMAL_ORDER = 5; - - /** Maximal aging triggering an attempt to balance the bracketing interval. */ - private static final int MAXIMAL_AGING = 2; - - /** Reduction factor for attempts to balance the bracketing interval. */ - private static final double REDUCTION_FACTOR = 1.0 / 16.0; - - /** Maximal order. */ - private final int maximalOrder; - - /** The kinds of solutions that the algorithm may accept. */ - private AllowedSolution allowed; - - /** - * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively) - */ - public BracketingNthOrderBrentSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER); - } - - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - * @param maximalOrder maximal order. - * @exception NumberIsTooSmallException if maximal order is lower than 2 - */ - public BracketingNthOrderBrentSolver(final double absoluteAccuracy, - final int maximalOrder) - throws NumberIsTooSmallException { - super(absoluteAccuracy); - if (maximalOrder < 2) { - throw new NumberIsTooSmallException(maximalOrder, 2, true); - } - this.maximalOrder = maximalOrder; - this.allowed = AllowedSolution.ANY_SIDE; - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param maximalOrder maximal order. - * @exception NumberIsTooSmallException if maximal order is lower than 2 - */ - public BracketingNthOrderBrentSolver(final double relativeAccuracy, - final double absoluteAccuracy, - final int maximalOrder) - throws NumberIsTooSmallException { - super(relativeAccuracy, absoluteAccuracy); - if (maximalOrder < 2) { - throw new NumberIsTooSmallException(maximalOrder, 2, true); - } - this.maximalOrder = maximalOrder; - this.allowed = AllowedSolution.ANY_SIDE; - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Function value accuracy. - * @param maximalOrder maximal order. - * @exception NumberIsTooSmallException if maximal order is lower than 2 - */ - public BracketingNthOrderBrentSolver(final double relativeAccuracy, - final double absoluteAccuracy, - final double functionValueAccuracy, - final int maximalOrder) - throws NumberIsTooSmallException { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); - if (maximalOrder < 2) { - throw new NumberIsTooSmallException(maximalOrder, 2, true); - } - this.maximalOrder = maximalOrder; - this.allowed = AllowedSolution.ANY_SIDE; - } - - /** Get the maximal order. - * @return maximal order - */ - public int getMaximalOrder() { - return maximalOrder; - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - // prepare arrays with the first points - final double[] x = new double[maximalOrder + 1]; - final double[] y = new double[maximalOrder + 1]; - x[0] = getMin(); - x[1] = getStartValue(); - x[2] = getMax(); - verifySequence(x[0], x[1], x[2]); - - // evaluate initial guess - y[1] = computeObjectiveValue(x[1]); - if (Precision.equals(y[1], 0.0, 1)) { - // return the initial guess if it is a perfect root. - return x[1]; - } - - // evaluate first endpoint - y[0] = computeObjectiveValue(x[0]); - if (Precision.equals(y[0], 0.0, 1)) { - // return the first endpoint if it is a perfect root. - return x[0]; - } - - int nbPoints; - int signChangeIndex; - if (y[0] * y[1] < 0) { - - // reduce interval if it brackets the root - nbPoints = 2; - signChangeIndex = 1; - - } else { - - // evaluate second endpoint - y[2] = computeObjectiveValue(x[2]); - if (Precision.equals(y[2], 0.0, 1)) { - // return the second endpoint if it is a perfect root. - return x[2]; - } - - if (y[1] * y[2] < 0) { - // use all computed point as a start sampling array for solving - nbPoints = 3; - signChangeIndex = 2; - } else { - throw new NoBracketingException(x[0], x[2], y[0], y[2]); - } - - } - - // prepare a work array for inverse polynomial interpolation - final double[] tmpX = new double[x.length]; - - // current tightest bracketing of the root - double xA = x[signChangeIndex - 1]; - double yA = y[signChangeIndex - 1]; - double absYA = FastMath.abs(yA); - int agingA = 0; - double xB = x[signChangeIndex]; - double yB = y[signChangeIndex]; - double absYB = FastMath.abs(yB); - int agingB = 0; - - // search loop - while (true) { - - // check convergence of bracketing interval - final double xTol = getAbsoluteAccuracy() + - getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB)); - if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) { - switch (allowed) { - case ANY_SIDE : - return absYA < absYB ? xA : xB; - case LEFT_SIDE : - return xA; - case RIGHT_SIDE : - return xB; - case BELOW_SIDE : - return (yA <= 0) ? xA : xB; - case ABOVE_SIDE : - return (yA < 0) ? xB : xA; - default : - // this should never happen - throw new MathInternalError(); - } - } - - // target for the next evaluation point - double targetY; - if (agingA >= MAXIMAL_AGING) { - // we keep updating the high bracket, try to compensate this - final int p = agingA - MAXIMAL_AGING; - final double weightA = (1 << p) - 1; - final double weightB = p + 1; - targetY = (weightA * yA - weightB * REDUCTION_FACTOR * yB) / (weightA + weightB); - } else if (agingB >= MAXIMAL_AGING) { - // we keep updating the low bracket, try to compensate this - final int p = agingB - MAXIMAL_AGING; - final double weightA = p + 1; - final double weightB = (1 << p) - 1; - targetY = (weightB * yB - weightA * REDUCTION_FACTOR * yA) / (weightA + weightB); - } else { - // bracketing is balanced, try to find the root itself - targetY = 0; - } - - // make a few attempts to guess a root, - double nextX; - int start = 0; - int end = nbPoints; - do { - - // guess a value for current target, using inverse polynomial interpolation - System.arraycopy(x, start, tmpX, start, end - start); - nextX = guessX(targetY, tmpX, y, start, end); - - if (!((nextX > xA) && (nextX < xB))) { - // the guessed root is not strictly inside of the tightest bracketing interval - - // the guessed root is either not strictly inside the interval or it - // is a NaN (which occurs when some sampling points share the same y) - // we try again with a lower interpolation order - if (signChangeIndex - start >= end - signChangeIndex) { - // we have more points before the sign change, drop the lowest point - ++start; - } else { - // we have more points after sign change, drop the highest point - --end; - } - - // we need to do one more attempt - nextX = Double.NaN; - - } - - } while (Double.isNaN(nextX) && (end - start > 1)); - - if (Double.isNaN(nextX)) { - // fall back to bisection - nextX = xA + 0.5 * (xB - xA); - start = signChangeIndex - 1; - end = signChangeIndex; - } - - // evaluate the function at the guessed root - final double nextY = computeObjectiveValue(nextX); - if (Precision.equals(nextY, 0.0, 1)) { - // we have found an exact root, since it is not an approximation - // we don't need to bother about the allowed solutions setting - return nextX; - } - - if ((nbPoints > 2) && (end - start != nbPoints)) { - - // we have been forced to ignore some points to keep bracketing, - // they are probably too far from the root, drop them from now on - nbPoints = end - start; - System.arraycopy(x, start, x, 0, nbPoints); - System.arraycopy(y, start, y, 0, nbPoints); - signChangeIndex -= start; - - } else if (nbPoints == x.length) { - - // we have to drop one point in order to insert the new one - nbPoints--; - - // keep the tightest bracketing interval as centered as possible - if (signChangeIndex >= (x.length + 1) / 2) { - // we drop the lowest point, we have to shift the arrays and the index - System.arraycopy(x, 1, x, 0, nbPoints); - System.arraycopy(y, 1, y, 0, nbPoints); - --signChangeIndex; - } - - } - - // insert the last computed point - //(by construction, we know it lies inside the tightest bracketing interval) - System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex); - x[signChangeIndex] = nextX; - System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex); - y[signChangeIndex] = nextY; - ++nbPoints; - - // update the bracketing interval - if (nextY * yA <= 0) { - // the sign change occurs before the inserted point - xB = nextX; - yB = nextY; - absYB = FastMath.abs(yB); - ++agingA; - agingB = 0; - } else { - // the sign change occurs after the inserted point - xA = nextX; - yA = nextY; - absYA = FastMath.abs(yA); - agingA = 0; - ++agingB; - - // update the sign change index - signChangeIndex++; - - } - - } - - } - - /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation. - * <p> - * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q - * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>), - * Q(y<sub>i</sub>) = x<sub>i</sub>. - * </p> - * @param targetY target value for y - * @param x reference points abscissas for interpolation, - * note that this array <em>is</em> modified during computation - * @param y reference points ordinates for interpolation - * @param start start index of the points to consider (inclusive) - * @param end end index of the points to consider (exclusive) - * @return guessed root (will be a NaN if two points share the same y) - */ - private double guessX(final double targetY, final double[] x, final double[] y, - final int start, final int end) { - - // compute Q Newton coefficients by divided differences - for (int i = start; i < end - 1; ++i) { - final int delta = i + 1 - start; - for (int j = end - 1; j > i; --j) { - x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]); - } - } - - // evaluate Q(targetY) - double x0 = 0; - for (int j = end - 1; j >= start; --j) { - x0 = x[j] + x0 * (targetY - y[j]); - } - - return x0; - - } - - /** {@inheritDoc} */ - public double solve(int maxEval, UnivariateFunction f, double min, - double max, AllowedSolution allowedSolution) - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - this.allowed = allowedSolution; - return super.solve(maxEval, f, min, max); - } - - /** {@inheritDoc} */ - public double solve(int maxEval, UnivariateFunction f, double min, - double max, double startValue, - AllowedSolution allowedSolution) - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - this.allowed = allowedSolution; - return super.solve(maxEval, f, min, max, startValue); - } - -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java deleted file mode 100644 index 0cc8750..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java +++ /dev/null @@ -1,244 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - - -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.Precision; - -/** - * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html";> - * Brent algorithm</a> for finding zeros of real univariate functions. - * The function should be continuous but not necessarily smooth. - * The {@code solve} method returns a zero {@code x} of the function {@code f} - * in the given interval {@code [a, b]} to within a tolerance - * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and - * {@code t} is the absolute accuracy. - * The given interval must bracket the root. - * <p> - * The reference implementation is given in chapter 4 of - * <blockquote> - * <b>Algorithms for Minimization Without Derivatives</b><br> - * <em>Richard P. Brent</em><br> - * Dover, 2002<br> - * </blockquote> - * </p> - * - * @see BaseAbstractUnivariateSolver - */ -public class BrentSolver extends AbstractUnivariateSolver { - - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver with default absolute accuracy (1e-6). - */ - public BrentSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public BrentSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public BrentSolver(double relativeAccuracy, - double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Function value accuracy. - * - * @see BaseAbstractUnivariateSolver#BaseAbstractUnivariateSolver(double,double,double) - */ - public BrentSolver(double relativeAccuracy, - double absoluteAccuracy, - double functionValueAccuracy) { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws NoBracketingException, - TooManyEvaluationsException, - NumberIsTooLargeException { - double min = getMin(); - double max = getMax(); - final double initial = getStartValue(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - - verifySequence(min, initial, max); - - // Return the initial guess if it is good enough. - double yInitial = computeObjectiveValue(initial); - if (FastMath.abs(yInitial) <= functionValueAccuracy) { - return initial; - } - - // Return the first endpoint if it is good enough. - double yMin = computeObjectiveValue(min); - if (FastMath.abs(yMin) <= functionValueAccuracy) { - return min; - } - - // Reduce interval if min and initial bracket the root. - if (yInitial * yMin < 0) { - return brent(min, initial, yMin, yInitial); - } - - // Return the second endpoint if it is good enough. - double yMax = computeObjectiveValue(max); - if (FastMath.abs(yMax) <= functionValueAccuracy) { - return max; - } - - // Reduce interval if initial and max bracket the root. - if (yInitial * yMax < 0) { - return brent(initial, max, yInitial, yMax); - } - - throw new NoBracketingException(min, max, yMin, yMax); - } - - /** - * Search for a zero inside the provided interval. - * This implementation is based on the algorithm described at page 58 of - * the book - * <blockquote> - * <b>Algorithms for Minimization Without Derivatives</b> - * <it>Richard P. Brent</it> - * Dover 0-486-41998-3 - * </blockquote> - * - * @param lo Lower bound of the search interval. - * @param hi Higher bound of the search interval. - * @param fLo Function value at the lower bound of the search interval. - * @param fHi Function value at the higher bound of the search interval. - * @return the value where the function is zero. - */ - private double brent(double lo, double hi, - double fLo, double fHi) { - double a = lo; - double fa = fLo; - double b = hi; - double fb = fHi; - double c = a; - double fc = fa; - double d = b - a; - double e = d; - - final double t = getAbsoluteAccuracy(); - final double eps = getRelativeAccuracy(); - - while (true) { - if (FastMath.abs(fc) < FastMath.abs(fb)) { - a = b; - b = c; - c = a; - fa = fb; - fb = fc; - fc = fa; - } - - final double tol = 2 * eps * FastMath.abs(b) + t; - final double m = 0.5 * (c - b); - - if (FastMath.abs(m) <= tol || - Precision.equals(fb, 0)) { - return b; - } - if (FastMath.abs(e) < tol || - FastMath.abs(fa) <= FastMath.abs(fb)) { - // Force bisection. - d = m; - e = d; - } else { - double s = fb / fa; - double p; - double q; - // The equality test (a == c) is intentional, - // it is part of the original Brent's method and - // it should NOT be replaced by proximity test. - if (a == c) { - // Linear interpolation. - p = 2 * m * s; - q = 1 - s; - } else { - // Inverse quadratic interpolation. - q = fa / fc; - final double r = fb / fc; - p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); - q = (q - 1) * (r - 1) * (s - 1); - } - if (p > 0) { - q = -q; - } else { - p = -p; - } - s = e; - e = d; - if (p >= 1.5 * m * q - FastMath.abs(tol * q) || - p >= FastMath.abs(0.5 * s * q)) { - // Inverse quadratic interpolation gives a value - // in the wrong direction, or progress is slow. - // Fall back to bisection. - d = m; - e = d; - } else { - d = p / q; - } - } - a = b; - fa = fb; - - if (FastMath.abs(d) > tol) { - b += d; - } else if (m > 0) { - b += tol; - } else { - b -= tol; - } - fb = computeObjectiveValue(b); - if ((fb > 0 && fc > 0) || - (fb <= 0 && fc <= 0)) { - c = a; - fc = fa; - d = b - a; - e = d; - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java deleted file mode 100644 index b9ae158..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java +++ /dev/null @@ -1,30 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; - - -/** - * Interface for (univariate real) rootfinding algorithms. - * Implementations will search for only one zero in the given interval. - * - * @deprecated as of 3.1, replaced by {@link UnivariateDifferentiableSolver} - */ -@Deprecated -public interface DifferentiableUnivariateSolver - extends BaseUnivariateSolver<DifferentiableUnivariateFunction> {} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java deleted file mode 100644 index bd3bc71..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java +++ /dev/null @@ -1,82 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.analysis.solvers; - - -/** - * Implements the <em>Illinois</em> method for root-finding (approximating - * a zero of a univariate real function). It is a modified - * {@link RegulaFalsiSolver <em>Regula Falsi</em>} method. - * - * <p>Like the <em>Regula Falsi</em> method, convergence is guaranteed by - * maintaining a bracketed solution. The <em>Illinois</em> method however, - * should converge much faster than the original <em>Regula Falsi</em> - * method. Furthermore, this implementation of the <em>Illinois</em> method - * should not suffer from the same implementation issues as the <em>Regula - * Falsi</em> method, which may fail to convergence in certain cases.</p> - * - * <p>The <em>Illinois</em> method assumes that the function is continuous, - * but not necessarily smooth.</p> - * - * <p>Implementation based on the following article: M. Dowell and P. Jarratt, - * <em>A modified regula falsi method for computing the root of an - * equation</em>, BIT Numerical Mathematics, volume 11, number 2, - * pages 168-174, Springer, 1971.</p> - * - * @since 3.0 - */ -public class IllinoisSolver extends BaseSecantSolver { - - /** Construct a solver with default accuracy (1e-6). */ - public IllinoisSolver() { - super(DEFAULT_ABSOLUTE_ACCURACY, Method.ILLINOIS); - } - - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public IllinoisSolver(final double absoluteAccuracy) { - super(absoluteAccuracy, Method.ILLINOIS); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public IllinoisSolver(final double relativeAccuracy, - final double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy, Method.ILLINOIS); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Maximum function value error. - */ - public IllinoisSolver(final double relativeAccuracy, - final double absoluteAccuracy, - final double functionValueAccuracy) { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.PEGASUS); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java deleted file mode 100644 index c127b42..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java +++ /dev/null @@ -1,390 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.complex.Complex; -import org.apache.commons.math3.complex.ComplexUtils; -import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.FastMath; - -/** - * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html";> - * Laguerre's Method</a> for root finding of real coefficient polynomials. - * For reference, see - * <blockquote> - * <b>A First Course in Numerical Analysis</b><br> - * ISBN 048641454X, chapter 8.<br> - * </blockquote> - * Laguerre's method is global in the sense that it can start with any initial - * approximation and be able to solve all roots from that point. - * The algorithm requires a bracketing condition. - * - * @since 1.2 - */ -public class LaguerreSolver extends AbstractPolynomialSolver { - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - /** Complex solver. */ - private final ComplexSolver complexSolver = new ComplexSolver(); - - /** - * Construct a solver with default accuracy (1e-6). - */ - public LaguerreSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public LaguerreSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public LaguerreSolver(double relativeAccuracy, - double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Function value accuracy. - */ - public LaguerreSolver(double relativeAccuracy, - double absoluteAccuracy, - double functionValueAccuracy) { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); - } - - /** - * {@inheritDoc} - */ - @Override - public double doSolve() - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - final double min = getMin(); - final double max = getMax(); - final double initial = getStartValue(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - - verifySequence(min, initial, max); - - // Return the initial guess if it is good enough. - final double yInitial = computeObjectiveValue(initial); - if (FastMath.abs(yInitial) <= functionValueAccuracy) { - return initial; - } - - // Return the first endpoint if it is good enough. - final double yMin = computeObjectiveValue(min); - if (FastMath.abs(yMin) <= functionValueAccuracy) { - return min; - } - - // Reduce interval if min and initial bracket the root. - if (yInitial * yMin < 0) { - return laguerre(min, initial, yMin, yInitial); - } - - // Return the second endpoint if it is good enough. - final double yMax = computeObjectiveValue(max); - if (FastMath.abs(yMax) <= functionValueAccuracy) { - return max; - } - - // Reduce interval if initial and max bracket the root. - if (yInitial * yMax < 0) { - return laguerre(initial, max, yInitial, yMax); - } - - throw new NoBracketingException(min, max, yMin, yMax); - } - - /** - * Find a real root in the given interval. - * - * Despite the bracketing condition, the root returned by - * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may - * not be a real zero inside {@code [min, max]}. - * For example, <code>p(x) = x<sup>3</sup> + 1,</code> - * with {@code min = -2}, {@code max = 2}, {@code initial = 0}. - * When it occurs, this code calls - * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)} - * in order to obtain all roots and picks up one real root. - * - * @param lo Lower bound of the search interval. - * @param hi Higher bound of the search interval. - * @param fLo Function value at the lower bound of the search interval. - * @param fHi Function value at the higher bound of the search interval. - * @return the point at which the function value is zero. - * @deprecated This method should not be part of the public API: It will - * be made private in version 4.0. - */ - @Deprecated - public double laguerre(double lo, double hi, - double fLo, double fHi) { - final Complex c[] = ComplexUtils.convertToComplex(getCoefficients()); - - final Complex initial = new Complex(0.5 * (lo + hi), 0); - final Complex z = complexSolver.solve(c, initial); - if (complexSolver.isRoot(lo, hi, z)) { - return z.getReal(); - } else { - double r = Double.NaN; - // Solve all roots and select the one we are seeking. - Complex[] root = complexSolver.solveAll(c, initial); - for (int i = 0; i < root.length; i++) { - if (complexSolver.isRoot(lo, hi, root[i])) { - r = root[i].getReal(); - break; - } - } - return r; - } - } - - /** - * Find all complex roots for the polynomial with the given - * coefficients, starting from the given initial value. - * <br/> - * Note: This method is not part of the API of {@link BaseUnivariateSolver}. - * - * @param coefficients Polynomial coefficients. - * @param initial Start value. - * @return the point at which the function value is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum number of evaluations is exceeded. - * @throws NullArgumentException if the {@code coefficients} is - * {@code null}. - * @throws NoDataException if the {@code coefficients} array is empty. - * @since 3.1 - */ - public Complex[] solveAllComplex(double[] coefficients, - double initial) - throws NullArgumentException, - NoDataException, - TooManyEvaluationsException { - setup(Integer.MAX_VALUE, - new PolynomialFunction(coefficients), - Double.NEGATIVE_INFINITY, - Double.POSITIVE_INFINITY, - initial); - return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients), - new Complex(initial, 0d)); - } - - /** - * Find a complex root for the polynomial with the given coefficients, - * starting from the given initial value. - * <br/> - * Note: This method is not part of the API of {@link BaseUnivariateSolver}. - * - * @param coefficients Polynomial coefficients. - * @param initial Start value. - * @return the point at which the function value is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum number of evaluations is exceeded. - * @throws NullArgumentException if the {@code coefficients} is - * {@code null}. - * @throws NoDataException if the {@code coefficients} array is empty. - * @since 3.1 - */ - public Complex solveComplex(double[] coefficients, - double initial) - throws NullArgumentException, - NoDataException, - TooManyEvaluationsException { - setup(Integer.MAX_VALUE, - new PolynomialFunction(coefficients), - Double.NEGATIVE_INFINITY, - Double.POSITIVE_INFINITY, - initial); - return complexSolver.solve(ComplexUtils.convertToComplex(coefficients), - new Complex(initial, 0d)); - } - - /** - * Class for searching all (complex) roots. - */ - private class ComplexSolver { - /** - * Check whether the given complex root is actually a real zero - * in the given interval, within the solver tolerance level. - * - * @param min Lower bound for the interval. - * @param max Upper bound for the interval. - * @param z Complex root. - * @return {@code true} if z is a real zero. - */ - public boolean isRoot(double min, double max, Complex z) { - if (isSequence(min, z.getReal(), max)) { - double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy()); - return (FastMath.abs(z.getImaginary()) <= tolerance) || - (z.abs() <= getFunctionValueAccuracy()); - } - return false; - } - - /** - * Find all complex roots for the polynomial with the given - * coefficients, starting from the given initial value. - * - * @param coefficients Polynomial coefficients. - * @param initial Start value. - * @return the point at which the function value is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum number of evaluations is exceeded. - * @throws NullArgumentException if the {@code coefficients} is - * {@code null}. - * @throws NoDataException if the {@code coefficients} array is empty. - */ - public Complex[] solveAll(Complex coefficients[], Complex initial) - throws NullArgumentException, - NoDataException, - TooManyEvaluationsException { - if (coefficients == null) { - throw new NullArgumentException(); - } - final int n = coefficients.length - 1; - if (n == 0) { - throw new NoDataException(LocalizedFormats.POLYNOMIAL); - } - // Coefficients for deflated polynomial. - final Complex c[] = new Complex[n + 1]; - for (int i = 0; i <= n; i++) { - c[i] = coefficients[i]; - } - - // Solve individual roots successively. - final Complex root[] = new Complex[n]; - for (int i = 0; i < n; i++) { - final Complex subarray[] = new Complex[n - i + 1]; - System.arraycopy(c, 0, subarray, 0, subarray.length); - root[i] = solve(subarray, initial); - // Polynomial deflation using synthetic division. - Complex newc = c[n - i]; - Complex oldc = null; - for (int j = n - i - 1; j >= 0; j--) { - oldc = c[j]; - c[j] = newc; - newc = oldc.add(newc.multiply(root[i])); - } - } - - return root; - } - - /** - * Find a complex root for the polynomial with the given coefficients, - * starting from the given initial value. - * - * @param coefficients Polynomial coefficients. - * @param initial Start value. - * @return the point at which the function value is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum number of evaluations is exceeded. - * @throws NullArgumentException if the {@code coefficients} is - * {@code null}. - * @throws NoDataException if the {@code coefficients} array is empty. - */ - public Complex solve(Complex coefficients[], Complex initial) - throws NullArgumentException, - NoDataException, - TooManyEvaluationsException { - if (coefficients == null) { - throw new NullArgumentException(); - } - - final int n = coefficients.length - 1; - if (n == 0) { - throw new NoDataException(LocalizedFormats.POLYNOMIAL); - } - - final double absoluteAccuracy = getAbsoluteAccuracy(); - final double relativeAccuracy = getRelativeAccuracy(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - - final Complex nC = new Complex(n, 0); - final Complex n1C = new Complex(n - 1, 0); - - Complex z = initial; - Complex oldz = new Complex(Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY); - while (true) { - // Compute pv (polynomial value), dv (derivative value), and - // d2v (second derivative value) simultaneously. - Complex pv = coefficients[n]; - Complex dv = Complex.ZERO; - Complex d2v = Complex.ZERO; - for (int j = n-1; j >= 0; j--) { - d2v = dv.add(z.multiply(d2v)); - dv = pv.add(z.multiply(dv)); - pv = coefficients[j].add(z.multiply(pv)); - } - d2v = d2v.multiply(new Complex(2.0, 0.0)); - - // Check for convergence. - final double tolerance = FastMath.max(relativeAccuracy * z.abs(), - absoluteAccuracy); - if ((z.subtract(oldz)).abs() <= tolerance) { - return z; - } - if (pv.abs() <= functionValueAccuracy) { - return z; - } - - // Now pv != 0, calculate the new approximation. - final Complex G = dv.divide(pv); - final Complex G2 = G.multiply(G); - final Complex H = G2.subtract(d2v.divide(pv)); - final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2)); - // Choose a denominator larger in magnitude. - final Complex deltaSqrt = delta.sqrt(); - final Complex dplus = G.add(deltaSqrt); - final Complex dminus = G.subtract(deltaSqrt); - final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus; - // Perturb z if denominator is zero, for instance, - // p(x) = x^3 + 1, z = 0. - if (denominator.equals(new Complex(0.0, 0.0))) { - z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy)); - oldz = new Complex(Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY); - } else { - oldz = z; - z = z.subtract(nC.divide(denominator)); - } - incrementEvaluationCount(); - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java deleted file mode 100644 index 06a7b6b..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java +++ /dev/null @@ -1,202 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; - -/** - * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html";> - * Muller's Method</a> for root finding of real univariate functions. For - * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, - * chapter 3. - * <p> - * Muller's method applies to both real and complex functions, but here we - * restrict ourselves to real functions. - * This class differs from {@link MullerSolver} in the way it avoids complex - * operations.</p> - * Muller's original method would have function evaluation at complex point. - * Since our f(x) is real, we have to find ways to avoid that. Bracketing - * condition is one way to go: by requiring bracketing in every iteration, - * the newly computed approximation is guaranteed to be real.</p> - * <p> - * Normally Muller's method converges quadratically in the vicinity of a - * zero, however it may be very slow in regions far away from zeros. For - * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use - * bisection as a safety backup if it performs very poorly.</p> - * <p> - * The formulas here use divided differences directly.</p> - * - * @since 1.2 - * @see MullerSolver2 - */ -public class MullerSolver extends AbstractUnivariateSolver { - - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver with default accuracy (1e-6). - */ - public MullerSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public MullerSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public MullerSolver(double relativeAccuracy, - double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - final double min = getMin(); - final double max = getMax(); - final double initial = getStartValue(); - - final double functionValueAccuracy = getFunctionValueAccuracy(); - - verifySequence(min, initial, max); - - // check for zeros before verifying bracketing - final double fMin = computeObjectiveValue(min); - if (FastMath.abs(fMin) < functionValueAccuracy) { - return min; - } - final double fMax = computeObjectiveValue(max); - if (FastMath.abs(fMax) < functionValueAccuracy) { - return max; - } - final double fInitial = computeObjectiveValue(initial); - if (FastMath.abs(fInitial) < functionValueAccuracy) { - return initial; - } - - verifyBracketing(min, max); - - if (isBracketing(min, initial)) { - return solve(min, initial, fMin, fInitial); - } else { - return solve(initial, max, fInitial, fMax); - } - } - - /** - * Find a real root in the given interval. - * - * @param min Lower bound for the interval. - * @param max Upper bound for the interval. - * @param fMin function value at the lower bound. - * @param fMax function value at the upper bound. - * @return the point at which the function value is zero. - * @throws TooManyEvaluationsException if the allowed number of calls to - * the function to be solved has been exhausted. - */ - private double solve(double min, double max, - double fMin, double fMax) - throws TooManyEvaluationsException { - final double relativeAccuracy = getRelativeAccuracy(); - final double absoluteAccuracy = getAbsoluteAccuracy(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - - // [x0, x2] is the bracketing interval in each iteration - // x1 is the last approximation and an interpolation point in (x0, x2) - // x is the new root approximation and new x1 for next round - // d01, d12, d012 are divided differences - - double x0 = min; - double y0 = fMin; - double x2 = max; - double y2 = fMax; - double x1 = 0.5 * (x0 + x2); - double y1 = computeObjectiveValue(x1); - - double oldx = Double.POSITIVE_INFINITY; - while (true) { - // Muller's method employs quadratic interpolation through - // x0, x1, x2 and x is the zero of the interpolating parabola. - // Due to bracketing condition, this parabola must have two - // real roots and we choose one in [x0, x2] to be x. - final double d01 = (y1 - y0) / (x1 - x0); - final double d12 = (y2 - y1) / (x2 - x1); - final double d012 = (d12 - d01) / (x2 - x0); - final double c1 = d01 + (x1 - x0) * d012; - final double delta = c1 * c1 - 4 * y1 * d012; - final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta)); - final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta)); - // xplus and xminus are two roots of parabola and at least - // one of them should lie in (x0, x2) - final double x = isSequence(x0, xplus, x2) ? xplus : xminus; - final double y = computeObjectiveValue(x); - - // check for convergence - final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); - if (FastMath.abs(x - oldx) <= tolerance || - FastMath.abs(y) <= functionValueAccuracy) { - return x; - } - - // Bisect if convergence is too slow. Bisection would waste - // our calculation of x, hopefully it won't happen often. - // the real number equality test x == x1 is intentional and - // completes the proximity tests above it - boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || - (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || - (x == x1); - // prepare the new bracketing interval for next iteration - if (!bisect) { - x0 = x < x1 ? x0 : x1; - y0 = x < x1 ? y0 : y1; - x2 = x > x1 ? x2 : x1; - y2 = x > x1 ? y2 : y1; - x1 = x; y1 = y; - oldx = x; - } else { - double xm = 0.5 * (x0 + x2); - double ym = computeObjectiveValue(xm); - if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) { - x2 = xm; y2 = ym; - } else { - x0 = xm; y0 = ym; - } - x1 = 0.5 * (x0 + x2); - y1 = computeObjectiveValue(x1); - oldx = Double.POSITIVE_INFINITY; - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java deleted file mode 100644 index 2a79c11..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java +++ /dev/null @@ -1,168 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; -import org.apache.commons.math3.util.FastMath; - -/** - * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html";> - * Muller's Method</a> for root finding of real univariate functions. For - * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, - * chapter 3. - * <p> - * Muller's method applies to both real and complex functions, but here we - * restrict ourselves to real functions. - * This class differs from {@link MullerSolver} in the way it avoids complex - * operations.</p> - * Except for the initial [min, max], it does not require bracketing - * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex - * number arises in the computation, we simply use its modulus as real - * approximation.</p> - * <p> - * Because the interval may not be bracketing, bisection alternative is - * not applicable here. However in practice our treatment usually works - * well, especially near real zeroes where the imaginary part of complex - * approximation is often negligible.</p> - * <p> - * The formulas here do not use divided differences directly.</p> - * - * @since 1.2 - * @see MullerSolver - */ -public class MullerSolver2 extends AbstractUnivariateSolver { - - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver with default accuracy (1e-6). - */ - public MullerSolver2() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public MullerSolver2(double absoluteAccuracy) { - super(absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public MullerSolver2(double relativeAccuracy, - double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException, - NumberIsTooLargeException, - NoBracketingException { - final double min = getMin(); - final double max = getMax(); - - verifyInterval(min, max); - - final double relativeAccuracy = getRelativeAccuracy(); - final double absoluteAccuracy = getAbsoluteAccuracy(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - - // x2 is the last root approximation - // x is the new approximation and new x2 for next round - // x0 < x1 < x2 does not hold here - - double x0 = min; - double y0 = computeObjectiveValue(x0); - if (FastMath.abs(y0) < functionValueAccuracy) { - return x0; - } - double x1 = max; - double y1 = computeObjectiveValue(x1); - if (FastMath.abs(y1) < functionValueAccuracy) { - return x1; - } - - if(y0 * y1 > 0) { - throw new NoBracketingException(x0, x1, y0, y1); - } - - double x2 = 0.5 * (x0 + x1); - double y2 = computeObjectiveValue(x2); - - double oldx = Double.POSITIVE_INFINITY; - while (true) { - // quadratic interpolation through x0, x1, x2 - final double q = (x2 - x1) / (x1 - x0); - final double a = q * (y2 - (1 + q) * y1 + q * y0); - final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; - final double c = (1 + q) * y2; - final double delta = b * b - 4 * a * c; - double x; - final double denominator; - if (delta >= 0.0) { - // choose a denominator larger in magnitude - double dplus = b + FastMath.sqrt(delta); - double dminus = b - FastMath.sqrt(delta); - denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus; - } else { - // take the modulus of (B +/- FastMath.sqrt(delta)) - denominator = FastMath.sqrt(b * b - delta); - } - if (denominator != 0) { - x = x2 - 2.0 * c * (x2 - x1) / denominator; - // perturb x if it exactly coincides with x1 or x2 - // the equality tests here are intentional - while (x == x1 || x == x2) { - x += absoluteAccuracy; - } - } else { - // extremely rare case, get a random number to skip it - x = min + FastMath.random() * (max - min); - oldx = Double.POSITIVE_INFINITY; - } - final double y = computeObjectiveValue(x); - - // check for convergence - final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); - if (FastMath.abs(x - oldx) <= tolerance || - FastMath.abs(y) <= functionValueAccuracy) { - return x; - } - - // prepare the next iteration - x0 = x1; - y0 = y1; - x1 = x2; - y1 = y2; - x2 = x; - y2 = y; - oldx = x; - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java deleted file mode 100644 index 4cf2688..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java +++ /dev/null @@ -1,92 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; -import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.exception.TooManyEvaluationsException; - -/** - * Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html";> - * Newton's Method</a> for finding zeros of real univariate differentiable - * functions. - * - * @since 3.1 - */ -public class NewtonRaphsonSolver extends AbstractUnivariateDifferentiableSolver { - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver. - */ - public NewtonRaphsonSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public NewtonRaphsonSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - - /** - * Find a zero near the midpoint of {@code min} and {@code max}. - * - * @param f Function to solve. - * @param min Lower bound for the interval. - * @param max Upper bound for the interval. - * @param maxEval Maximum number of evaluations. - * @return the value where the function is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum evaluation count is exceeded. - * @throws org.apache.commons.math3.exception.NumberIsTooLargeException - * if {@code min >= max}. - */ - @Override - public double solve(int maxEval, final UnivariateDifferentiableFunction f, - final double min, final double max) - throws TooManyEvaluationsException { - return super.solve(maxEval, f, UnivariateSolverUtils.midpoint(min, max)); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException { - final double startValue = getStartValue(); - final double absoluteAccuracy = getAbsoluteAccuracy(); - - double x0 = startValue; - double x1; - while (true) { - final DerivativeStructure y0 = computeObjectiveValueAndDerivative(x0); - x1 = x0 - (y0.getValue() / y0.getPartialDerivative(1)); - if (FastMath.abs(x1 - x0) <= absoluteAccuracy) { - return x1; - } - - x0 = x1; - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java deleted file mode 100644 index 3ba7bf2..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java +++ /dev/null @@ -1,92 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.exception.TooManyEvaluationsException; - -/** - * Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html";> - * Newton's Method</a> for finding zeros of real univariate functions. - * <p> - * The function should be continuous but not necessarily smooth.</p> - * - * @deprecated as of 3.1, replaced by {@link NewtonRaphsonSolver} - */ -@Deprecated -public class NewtonSolver extends AbstractDifferentiableUnivariateSolver { - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver. - */ - public NewtonSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public NewtonSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - - /** - * Find a zero near the midpoint of {@code min} and {@code max}. - * - * @param f Function to solve. - * @param min Lower bound for the interval. - * @param max Upper bound for the interval. - * @param maxEval Maximum number of evaluations. - * @return the value where the function is zero. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the maximum evaluation count is exceeded. - * @throws org.apache.commons.math3.exception.NumberIsTooLargeException - * if {@code min >= max}. - */ - @Override - public double solve(int maxEval, final DifferentiableUnivariateFunction f, - final double min, final double max) - throws TooManyEvaluationsException { - return super.solve(maxEval, f, UnivariateSolverUtils.midpoint(min, max)); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException { - final double startValue = getStartValue(); - final double absoluteAccuracy = getAbsoluteAccuracy(); - - double x0 = startValue; - double x1; - while (true) { - x1 = x0 - (computeObjectiveValue(x0) / computeDerivativeObjectiveValue(x0)); - if (FastMath.abs(x1 - x0) <= absoluteAccuracy) { - return x1; - } - - x0 = x1; - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java deleted file mode 100644 index 0d80895..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java +++ /dev/null @@ -1,84 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.analysis.solvers; - -/** - * Implements the <em>Pegasus</em> method for root-finding (approximating - * a zero of a univariate real function). It is a modified - * {@link RegulaFalsiSolver <em>Regula Falsi</em>} method. - * - * <p>Like the <em>Regula Falsi</em> method, convergence is guaranteed by - * maintaining a bracketed solution. The <em>Pegasus</em> method however, - * should converge much faster than the original <em>Regula Falsi</em> - * method. Furthermore, this implementation of the <em>Pegasus</em> method - * should not suffer from the same implementation issues as the <em>Regula - * Falsi</em> method, which may fail to convergence in certain cases. Also, - * the <em>Pegasus</em> method should converge faster than the - * {@link IllinoisSolver <em>Illinois</em>} method, another <em>Regula - * Falsi</em>-based method.</p> - * - * <p>The <em>Pegasus</em> method assumes that the function is continuous, - * but not necessarily smooth.</p> - * - * <p>Implementation based on the following article: M. Dowell and P. Jarratt, - * <em>The "Pegasus" method for computing the root of an equation</em>, - * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer, - * 1972.</p> - * - * @since 3.0 - */ -public class PegasusSolver extends BaseSecantSolver { - - /** Construct a solver with default accuracy (1e-6). */ - public PegasusSolver() { - super(DEFAULT_ABSOLUTE_ACCURACY, Method.PEGASUS); - } - - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public PegasusSolver(final double absoluteAccuracy) { - super(absoluteAccuracy, Method.PEGASUS); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public PegasusSolver(final double relativeAccuracy, - final double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy, Method.PEGASUS); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Maximum function value error. - */ - public PegasusSolver(final double relativeAccuracy, - final double absoluteAccuracy, - final double functionValueAccuracy) { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.PEGASUS); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java deleted file mode 100644 index c21f076..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java +++ /dev/null @@ -1,28 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; - -/** - * Interface for (polynomial) root-finding algorithms. - * Implementations will search for only one zero in the given interval. - * - * @since 3.0 - */ -public interface PolynomialSolver - extends BaseUnivariateSolver<PolynomialFunction> {} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java deleted file mode 100644 index cfb7055..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java +++ /dev/null @@ -1,94 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.analysis.solvers; - -/** - * Implements the <em>Regula Falsi</em> or <em>False position</em> method for - * root-finding (approximating a zero of a univariate real function). It is a - * modified {@link SecantSolver <em>Secant</em>} method. - * - * <p>The <em>Regula Falsi</em> method is included for completeness, for - * testing purposes, for educational purposes, for comparison to other - * algorithms, etc. It is however <strong>not</strong> intended to be used - * for actual problems, as one of the bounds often remains fixed, resulting - * in very slow convergence. Instead, one of the well-known modified - * <em>Regula Falsi</em> algorithms can be used ({@link IllinoisSolver - * <em>Illinois</em>} or {@link PegasusSolver <em>Pegasus</em>}). These two - * algorithms solve the fundamental issues of the original <em>Regula - * Falsi</em> algorithm, and greatly out-performs it for most, if not all, - * (practical) functions. - * - * <p>Unlike the <em>Secant</em> method, the <em>Regula Falsi</em> guarantees - * convergence, by maintaining a bracketed solution. Note however, that due to - * the finite/limited precision of Java's {@link Double double} type, which is - * used in this implementation, the algorithm may get stuck in a situation - * where it no longer makes any progress. Such cases are detected and result - * in a {@code ConvergenceException} exception being thrown. In other words, - * the algorithm theoretically guarantees convergence, but the implementation - * does not.</p> - * - * <p>The <em>Regula Falsi</em> method assumes that the function is continuous, - * but not necessarily smooth.</p> - * - * <p>Implementation based on the following article: M. Dowell and P. Jarratt, - * <em>A modified regula falsi method for computing the root of an - * equation</em>, BIT Numerical Mathematics, volume 11, number 2, - * pages 168-174, Springer, 1971.</p> - * - * @since 3.0 - */ -public class RegulaFalsiSolver extends BaseSecantSolver { - - /** Construct a solver with default accuracy (1e-6). */ - public RegulaFalsiSolver() { - super(DEFAULT_ABSOLUTE_ACCURACY, Method.REGULA_FALSI); - } - - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public RegulaFalsiSolver(final double absoluteAccuracy) { - super(absoluteAccuracy, Method.REGULA_FALSI); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public RegulaFalsiSolver(final double relativeAccuracy, - final double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy, Method.REGULA_FALSI); - } - - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - * @param functionValueAccuracy Maximum function value error. - */ - public RegulaFalsiSolver(final double relativeAccuracy, - final double absoluteAccuracy, - final double functionValueAccuracy) { - super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.REGULA_FALSI); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java deleted file mode 100644 index d83f595..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java +++ /dev/null @@ -1,142 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.solvers; - -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.exception.NoBracketingException; -import org.apache.commons.math3.exception.TooManyEvaluationsException; - -/** - * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html";> - * Ridders' Method</a> for root finding of real univariate functions. For - * reference, see C. Ridders, <i>A new algorithm for computing a single root - * of a real continuous function </i>, IEEE Transactions on Circuits and - * Systems, 26 (1979), 979 - 980. - * <p> - * The function should be continuous but not necessarily smooth.</p> - * - * @since 1.2 - */ -public class RiddersSolver extends AbstractUnivariateSolver { - /** Default absolute accuracy. */ - private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; - - /** - * Construct a solver with default accuracy (1e-6). - */ - public RiddersSolver() { - this(DEFAULT_ABSOLUTE_ACCURACY); - } - /** - * Construct a solver. - * - * @param absoluteAccuracy Absolute accuracy. - */ - public RiddersSolver(double absoluteAccuracy) { - super(absoluteAccuracy); - } - /** - * Construct a solver. - * - * @param relativeAccuracy Relative accuracy. - * @param absoluteAccuracy Absolute accuracy. - */ - public RiddersSolver(double relativeAccuracy, - double absoluteAccuracy) { - super(relativeAccuracy, absoluteAccuracy); - } - - /** - * {@inheritDoc} - */ - @Override - protected double doSolve() - throws TooManyEvaluationsException, - NoBracketingException { - double min = getMin(); - double max = getMax(); - // [x1, x2] is the bracketing interval in each iteration - // x3 is the midpoint of [x1, x2] - // x is the new root approximation and an endpoint of the new interval - double x1 = min; - double y1 = computeObjectiveValue(x1); - double x2 = max; - double y2 = computeObjectiveValue(x2); - - // check for zeros before verifying bracketing - if (y1 == 0) { - return min; - } - if (y2 == 0) { - return max; - } - verifyBracketing(min, max); - - final double absoluteAccuracy = getAbsoluteAccuracy(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - final double relativeAccuracy = getRelativeAccuracy(); - - double oldx = Double.POSITIVE_INFINITY; - while (true) { - // calculate the new root approximation - final double x3 = 0.5 * (x1 + x2); - final double y3 = computeObjectiveValue(x3); - if (FastMath.abs(y3) <= functionValueAccuracy) { - return x3; - } - final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing - final double correction = (FastMath.signum(y2) * FastMath.signum(y3)) * - (x3 - x1) / FastMath.sqrt(delta); - final double x = x3 - correction; // correction != 0 - final double y = computeObjectiveValue(x); - - // check for convergence - final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); - if (FastMath.abs(x - oldx) <= tolerance) { - return x; - } - if (FastMath.abs(y) <= functionValueAccuracy) { - return x; - } - - // prepare the new interval for next iteration - // Ridders' method guarantees x1 < x < x2 - if (correction > 0.0) { // x1 < x < x3 - if (FastMath.signum(y1) + FastMath.signum(y) == 0.0) { - x2 = x; - y2 = y; - } else { - x1 = x; - x2 = x3; - y1 = y; - y2 = y3; - } - } else { // x3 < x < x2 - if (FastMath.signum(y2) + FastMath.signum(y) == 0.0) { - x1 = x; - y1 = y; - } else { - x1 = x3; - x2 = x; - y1 = y3; - y2 = y; - } - } - oldx = x; - } - } -}
