diff --git a/src/java/org/apache/commons/math/analysis/BisectionSolver.java b/src/java/org/apache/commons/math/analysis/BisectionSolver.java new file mode 100644 index 000000000..2222b2ab5 --- /dev/null +++ b/src/java/org/apache/commons/math/analysis/BisectionSolver.java @@ -0,0 +1,139 @@ +/* ==================================================================== + * The Apache Software License, Version 1.1 + * + * Copyright (c) 2003 The Apache Software Foundation. All rights + * reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * 3. The end-user documentation included with the redistribution, if + * any, must include the following acknowlegement: + * "This product includes software developed by the + * Apache Software Foundation (http://www.apache.org/)." + * Alternately, this acknowlegement may appear in the software itself, + * if and wherever such third-party acknowlegements normally appear. + * + * 4. The names "The Jakarta Project", "Commons", and "Apache Software + * Foundation" must not be used to endorse or promote products derived + * from this software without prior written permission. For written + * permission, please contact apache@apache.org. + * + * 5. Products derived from this software may not be called "Apache" + * nor may "Apache" appear in their names without prior written + * permission of the Apache Software Foundation. + * + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF + * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT + * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * ==================================================================== + * + * This software consists of voluntary contributions made by many + * individuals on behalf of the Apache Software Foundation. For more + * information on the Apache Software Foundation, please see + * . + */ +package org.apache.commons.math.analysis; + +import org.apache.commons.math.MathException; + +/** + * Provide the bisection algorithm for solving for zeros of real univariate + * functions. It will only search for one zero in the given interval. The + * function is supposed to be continuous but not necessarily smooth. + * + * @author Brent Worden + */ +public class BisectionSolver extends UnivariateRealSolverImpl { + /** + * Construct a solver for the given function. + * @param f function to solve. + */ + public BisectionSolver(UnivariateRealFunction f) { + super(f, 100, 1E-6); + } + + /** + * Solve for a zero in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @param initial the start value to use (ignored). + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. + */ + public double solve(double min, double max, double initial) + throws MathException { + + return solve(min, max); + } + + /** + * Solve for a zero root in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. + */ + public double solve(double min, double max) throws MathException { + clearResult(); + + double m; + double fm; + double fmin; + + int i = 0; + while (i < maximalIterationCount) { + m = midpoint(min, max); + fmin = f.value(min); + fm = f.value(m); + + if (fm * fmin > 0.0) { + // max and m bracket the root. + min = m; + fmin = fm; + } else { + // min and m bracket the root. + max = m; + } + + if (Math.abs(max - min) <= absoluteAccuracy) { + m = midpoint(min, max); + setResult(m, i); + return m; + } + ++i; + } + + throw new MathException("Maximal iteration number exceeded"); + } + + /** + * Compute the midpoint of two values. + * @param a first value. + * @param b second value. + * @return the midpoint. + */ + public static double midpoint(double a, double b) { + return (a + b) * .5; + } +} diff --git a/src/java/org/apache/commons/math/analysis/BrentSolver.java b/src/java/org/apache/commons/math/analysis/BrentSolver.java index 6541c3e95..4941f300f 100644 --- a/src/java/org/apache/commons/math/analysis/BrentSolver.java +++ b/src/java/org/apache/commons/math/analysis/BrentSolver.java @@ -64,16 +64,36 @@ import org.apache.commons.math.MathException; * @author pietsch at apache.org */ public class BrentSolver extends UnivariateRealSolverImpl { - - private UnivariateRealFunction f; - + /** + * Construct a solver for the given function. + * @param f function to solve. + */ public BrentSolver(UnivariateRealFunction f) { - super(100, 1E-6); - this.f = f; + super(f, 100, 1E-6); } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#solve(double, double) + /** + * Solve for a zero in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @param initial the start value to use (ignored). + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. + */ + public double solve(double min, double max, double initial) + throws MathException { + + return solve(min, max); + } + + /** + * Solve for a zero root in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. */ public double solve(double min, double max) throws MathException { clearResult(); @@ -116,21 +136,8 @@ public class BrentSolver extends UnivariateRealSolverImpl { setResult(x1, i); return result; } -// System.out.println( -// " x0=" -// + x0 -// + " y0=" -// + y0 -// + " x1=" -// + x1 -// + " y1=" -// + y1 -// + " x2=" -// + x2+" y2="+y2); -// System.out.println(" dx="+dx+" delta: "+delta+" olddelta: "+oldDelta); if (Math.abs(oldDelta) < tolerance || Math.abs(y0) <= Math.abs(y1)) { -// System.out.println("bisection"); // Force bisection. delta = 0.5 * dx; oldDelta = delta; @@ -140,12 +147,10 @@ public class BrentSolver extends UnivariateRealSolverImpl { double p1; if (x0 == x2) { // Linear interpolation. -// System.out.println("linear"); p = dx * r3; p1 = 1.0 - r3; } else { // Inverse quadratic interpolation. -// System.out.println("invers quad"); double r1 = y0 / y2; double r2 = y1 / y2; p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0)); @@ -156,14 +161,11 @@ public class BrentSolver extends UnivariateRealSolverImpl { } else { p = -p; } -// System.out.println(" p="+p+" p1="+p1+" qq: "+(1.5 * dx * p1 - Math.abs(tolerance * p1))); -// System.out.println(" p="+p+" q: "+p1+" ad="+Math.abs(0.5 * oldDelta * p1)); if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) || p >= Math.abs(0.5 * oldDelta * p1)) { // Inverse quadratic interpolation gives a value // in the wrong direction, or progress is slow. // Fall back to bisection. -// System.out.println("bisection fallback"); delta = 0.5 * dx; oldDelta = delta; } else { @@ -178,9 +180,9 @@ public class BrentSolver extends UnivariateRealSolverImpl { if (Math.abs(delta) > tolerance) { x1 = x1 + delta; } else if (dx > 0.0) { - x1 = x1 + 0.5*tolerance; + x1 = x1 + 0.5 * tolerance; } else if (dx <= 0.0) { - x1 = x1 - 0.5*tolerance; + x1 = x1 - 0.5 * tolerance; } y1 = f.value(x1); if ((y1 > 0) == (y2 > 0)) { @@ -193,5 +195,4 @@ public class BrentSolver extends UnivariateRealSolverImpl { } throw new MathException("Maximal iteration number exceeded."); } - } diff --git a/src/java/org/apache/commons/math/analysis/RootFinding.java b/src/java/org/apache/commons/math/analysis/RootFinding.java index 0c7de3f7c..4ce8ccfc2 100644 --- a/src/java/org/apache/commons/math/analysis/RootFinding.java +++ b/src/java/org/apache/commons/math/analysis/RootFinding.java @@ -53,6 +53,8 @@ */ package org.apache.commons.math.analysis; +import org.apache.commons.math.MathException; + /** * Utility class comprised of root finding techniques. * @@ -77,14 +79,17 @@ public class RootFinding { * @param function the function * @param initial midpoint of the returned range. * @param lowerBound for numerical safety, a never is less than this value. - * @param upperBound for numerical safety, b never is greater than this value. + * @param upperBound for numerical safety, b never is greater than this + * value. * @return a two element array holding {a, b}. + * @throws MathException if a root can not be bracketted. */ - public static double[] bracket(UnivariateFunction function, + public static double[] bracket(UnivariateRealFunction function, double initial, double lowerBound, - double upperBound) { - return bracket( function, initial, lowerBound, upperBound, Integer.MAX_VALUE ) ; + double upperBound) throws MathException { + return bracket( function, initial, lowerBound, upperBound, + Integer.MAX_VALUE ) ; } /** @@ -95,15 +100,18 @@ public class RootFinding { * @param function the function * @param initial midpoint of the returned range. * @param lowerBound for numerical safety, a never is less than this value. - * @param upperBound for numerical safety, b never is greater than this value. - * @param maximumIterations to guard against infinite looping, maximum number of iterations to perform + * @param upperBound for numerical safety, b never is greater than this + * value. + * @param maximumIterations to guard against infinite looping, maximum + * number of iterations to perform * @return a two element array holding {a, b}. + * @throws MathException if a root can not be bracketted. */ - public static double[] bracket(UnivariateFunction function, + public static double[] bracket(UnivariateRealFunction function, double initial, double lowerBound, double upperBound, - int maximumIterations) { + int maximumIterations) throws MathException { double a = initial; double b = initial; double fa; @@ -113,52 +121,11 @@ public class RootFinding { do { a = Math.max(a - 1.0, lowerBound); b = Math.min(b + 1.0, upperBound); - fa = function.evaluate(a); - fb = function.evaluate(b); + fa = function.value(a); + fb = function.value(b); numIterations += 1 ; } while ( (fa * fb > 0.0) && ( numIterations < maximumIterations ) ); return new double[]{a, b}; } - - /** - * For a function, f, this method returns a root c that lies between a and - * b, and satisfies f(c) = 0. - * - * @param function the function - * @param a lower (or upper) bound of a root - * @param b upper (or lower) bound of a root - * @return a root of f - */ - public static double bisection(UnivariateFunction function, - double a, - double b) { - double m; - double fm; - double fa; - - if ( b < a ) { - double xtemp = a ; - a = b ; - b = xtemp ; - } - - fa = function.evaluate(a); - - while(Math.abs(a - b) > EPSILON) { - m = (a + b) * 0.5; // midpoint - fm = function.evaluate(m); - - if(fm * fa > 0.0) { - // b and m bracket the root. - a = m; - fa = fm; - } else { - // a and m bracket the root. - b = m; - } - } - - return (a + b) * 0.5; - } } diff --git a/src/java/org/apache/commons/math/analysis/SecantSolver.java b/src/java/org/apache/commons/math/analysis/SecantSolver.java index 140261bcf..cb7158112 100644 --- a/src/java/org/apache/commons/math/analysis/SecantSolver.java +++ b/src/java/org/apache/commons/math/analysis/SecantSolver.java @@ -66,16 +66,36 @@ import org.apache.commons.math.MathException; * @author pietsch at apache.org */ public class SecantSolver extends UnivariateRealSolverImpl { - - private UnivariateRealFunction f; - + /** + * Construct a solver for the given function. + * @param f function to solve. + */ public SecantSolver(UnivariateRealFunction f) { - super(100, 1E-6); - this.f = f; + super(f, 100, 1E-6); } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#solve(double, double) + /** + * Solve for a zero in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @param initial the start value to use (ignored). + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. + */ + public double solve(double min, double max, double initial) + throws MathException { + + return solve(min, max); + } + + /** + * Solve for a zero root in the given interval. + * @param min the lower bound for the interval. + * @param max the upper bound for the interval. + * @return the value where the function is zero + * @throws MathException if the iteration count was exceeded or the + * solver detects convergence problems otherwise. */ public double solve(double min, double max) throws MathException { clearResult(); @@ -88,7 +108,7 @@ public class SecantSolver extends UnivariateRealSolverImpl { double x1 = max; double y0 = f.value(x0); double y1 = f.value(x1); - if ((y0>0)== (y1>0)) { + if ((y0 > 0) == (y1 > 0)) { throw new MathException("Interval doesn't bracket a zero."); } double x2 = x0; @@ -117,21 +137,19 @@ public class SecantSolver extends UnivariateRealSolverImpl { if (Math.abs(y1) > Math.abs(y0)) { // Function value increased in last iteration. Force bisection. delta = 0.5 * oldDelta; -// System.out.println("Forced Bisection"); } else { delta = (x0 - x1) / (1 - y0 / y1); - // System.out.println("delta=" + delta + " olddelta=" + oldDelta); if (delta / oldDelta > 1) { - // New approximation falls outside bracket. Fall back to bisection. + // New approximation falls outside bracket. + // Fall back to bisection. delta = 0.5 * oldDelta; -// System.out.println("Fallback Bisection"); } } x0 = x1; y0 = y1; x1 = x1 + delta; y1 = f.value(x1); - if ((y1>0) == (y2>0)) { + if ((y1 > 0) == (y2 > 0)) { // New bracket is (x0,x1). x2 = x0; y2 = y0; diff --git a/src/java/org/apache/commons/math/analysis/UnivariateRealSolver.java b/src/java/org/apache/commons/math/analysis/UnivariateRealSolver.java index c8669d451..c03e8f03a 100644 --- a/src/java/org/apache/commons/math/analysis/UnivariateRealSolver.java +++ b/src/java/org/apache/commons/math/analysis/UnivariateRealSolver.java @@ -66,23 +66,27 @@ public interface UnivariateRealSolver { /** * Set the upper limit for the number of iterations. + * * Usually a high iteration count indicates convergence problems. However, * the "reasonable value" varies widely for different solvers, users are - * advised to use the default value supplied by the solver. + * advised to use the default value supplied by the solver. + * * An exception will be thrown if the number is exceeded. * - * @param count + * @param count maximum number of iterations */ public void setMaximalIterationCount(int count); /** * Get the upper limit for the number of iterations. + * * @return the actual upper limit */ public int getMaximalIterationCount(); /** * Reset the upper limit for the number of iterations to the default. + * * The default value is supplied by the solver implementation. * * @see #setMaximalIterationCount(int) @@ -91,41 +95,50 @@ public interface UnivariateRealSolver { /** * Set the absolute accuracy. + * * The default is usually choosen so taht roots in the interval - * -10..-0.1 and +0.1..+10 can be found wit a reasonable accuracy. If the expected - * absolute value of your roots is of much smaller magnitude, set this to a smaller - * value. - * Solvers are advised to do a plausibility check with the relative accuracy, but - * clients should not rely on this. + * -10..-0.1 and +0.1..+10 can be found wit a reasonable accuracy. If the + * expected absolute value of your roots is of much smaller magnitude, set + * this to a smaller value. + * + * Solvers are advised to do a plausibility check with the relative + * accuracy, but clients should not rely on this. + * * @param accuracy the accuracy. - * @throws MathException if the accuracy can't be achieved by the solver or is - * otherwise deemed unreasonable. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ public void setAbsoluteAccuracy(double accuracy) throws MathException; /** * Get the actual absolute accuracy. + * * @return the accuracy */ public double getAbsoluteAccuracy(); /** * Reset the absolute accuracy to the default. + * * The default value is provided by the solver implementation. */ public void resetAbsoluteAccuracy(); /** * Set the relative accuracy. - * This is used to stop iterations if the absolute accuracy can't be achieved - * due to large values or short mantissa length. - * If this should be the primary criterium for convergence rather then a safety - * measure, set the absolute accuracy to a ridiculously small value, like 1E-1000. + * + * This is used to stop iterations if the absolute accuracy can't be + * achieved due to large values or short mantissa length. + * + * If this should be the primary criterium for convergence rather then a + * safety measure, set the absolute accuracy to a ridiculously small value, + * like 1E-1000. + * * @param accuracy the relative accuracy. - * @throws MathException if the accuracy can't be achieved by the solver or is - * otherwise deemed unreasonable. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ - public void setRelativeAccuracy(double Accuracy) throws MathException; + public void setRelativeAccuracy(double accuracy) throws MathException; /** * Get the actual relative accuracy. @@ -141,14 +154,18 @@ public interface UnivariateRealSolver { /** * Set the function value accuracy. + * * This is used to determine whan an evaluated function value or some other * value which is used as divisor is zero. - * This is a safety guard and it shouldn't be necesary to change this in general. + * + * This is a safety guard and it shouldn't be necesary to change this in + * general. + * * @param accuracy the accuracy. - * @throws MathException if the accuracy can't be achieved by the solver or is - * otherwise deemed unreasonable. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ - public void setFunctionValueAccuracy(double Accuracy) throws MathException; + public void setFunctionValueAccuracy(double accuracy) throws MathException; /** * Get the actual function value accuracy. diff --git a/src/java/org/apache/commons/math/analysis/UnivariateRealSolverImpl.java b/src/java/org/apache/commons/math/analysis/UnivariateRealSolverImpl.java index 3c825b868..cc3ede12b 100644 --- a/src/java/org/apache/commons/math/analysis/UnivariateRealSolverImpl.java +++ b/src/java/org/apache/commons/math/analysis/UnivariateRealSolverImpl.java @@ -65,24 +65,57 @@ import org.apache.commons.math.MathException; public abstract class UnivariateRealSolverImpl implements UnivariateRealSolver { + /** Maximum absolute error. */ protected double absoluteAccuracy; + + /** Maximum relative error. */ protected double relativeAccuracy; + + /** Maximum error of function. */ protected double functionValueAccuracy; + + /** Maximum number of iterations. */ protected int maximalIterationCount; + /** Default maximum absolute error. */ protected double defaultAbsoluteAccuracy; + + /** Default maximum relative error. */ protected double defaultRelativeAccuracy; + + /** Default maximum error of function. */ protected double defaultFunctionValueAccuracy; + + /** Default maximum number of iterations. */ protected int defaultMaximalIterationCount; + /** Indicates where a root has been computed. */ protected boolean resultComputed = false; + + /** The last computed root. */ protected double result; + // Mainly for test framework. + /** The last iteration count. */ protected int iterationCount; + /** The function to solve. */ + protected UnivariateRealFunction f; + + /** + * Construct a solver with given iteration count and accuracy. + * @param f the function to solve. + * @param defaultAbsoluteAccuracy maximum absolue error. + * @param defaultMaximalIterationCount maximum number of iterations. + */ protected UnivariateRealSolverImpl( + UnivariateRealFunction f, int defaultMaximalIterationCount, double defaultAbsoluteAccuracy) { + + super(); + + this.f = f; this.defaultAbsoluteAccuracy = defaultAbsoluteAccuracy; this.defaultRelativeAccuracy = 1E-14; this.defaultFunctionValueAccuracy = 1E-15; @@ -93,48 +126,39 @@ public abstract class UnivariateRealSolverImpl this.maximalIterationCount = defaultMaximalIterationCount; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#solve(double, double) - */ - public double solve(double min, double max) throws MathException { - throw new UnsupportedOperationException(); - } - - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#solve(double, double, double) - */ - public double solve(double min, double max, double startValue) - throws MathException { - throw new UnsupportedOperationException(); - } - - /* - * Get result of last solver run. - * @see org.apache.commons.math.UnivariateRealSolver#getResult() + /** + * Access the last computed root. + * @return the last computed root. + * @throws MathException if no root has been computed. */ public double getResult() throws MathException { if (resultComputed) { return result; } else { + // TODO: could this be an IllegalStateException instead? throw new MathException("No result available"); } } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#getIterationCount() + /** + * Access the last iteration count. + * @return the last iteration count. + * @throws MathException if no root has been computed. + * */ public int getIterationCount() throws MathException { if (resultComputed) { return iterationCount; } else { + // TODO: could this be an IllegalStateException instead? throw new MathException("No result available"); } } - /* + /** * Convenience function for implementations. * @param result the result to set - * @param iteratinCount the iteration count to set + * @param iterationCount the iteration count to set */ protected final void setResult(double result, int iterationCount) { this.result = result; @@ -142,97 +166,116 @@ public abstract class UnivariateRealSolverImpl this.resultComputed = true; } - /* + /** * Convenience function for implementations. */ protected final void clearResult() { this.resultComputed = false; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#setAccuracy(double) + /** + * Set the absolute accuracy. + * + * @param accuracy the accuracy. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ public void setAbsoluteAccuracy(double accuracy) throws MathException { absoluteAccuracy = accuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#getAccuracy() + /** + * Get the actual absolute accuracy. + * + * @return the accuracy */ public double getAbsoluteAccuracy() { return absoluteAccuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#resetAbsoluteAccuracy() + /** + * Reset the absolute accuracy to the default. */ public void resetAbsoluteAccuracy() { absoluteAccuracy = defaultAbsoluteAccuracy; } - /* Set maximum iteration count. - * @see org.apache.commons.math.UnivariateRealSolver#setMaximalIterationCount(int) + /** + * Set the upper limit for the number of iterations. + * + * @param count maximum number of iterations */ public void setMaximalIterationCount(int count) { maximalIterationCount = count; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#getMaximalIterationCount() + /** + * Get the upper limit for the number of iterations. + * + * @return the actual upper limit */ public int getMaximalIterationCount() { return maximalIterationCount; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#resetMaximalIterationCount() + /** + * Reset the upper limit for the number of iterations to the default. */ public void resetMaximalIterationCount() { maximalIterationCount = defaultMaximalIterationCount; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#setRelativeAccuracy(double) + /** + * Set the relative accuracy. + * + * @param accuracy the relative accuracy. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ public void setRelativeAccuracy(double accuracy) throws MathException { relativeAccuracy = accuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#getRelativeAccuracy() + /** + * Get the actual relative accuracy. + * @return the accuracy */ public double getRelativeAccuracy() { return relativeAccuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#resetRelativeAccuracy() + /** + * Reset the relative accuracy to the default. */ public void resetRelativeAccuracy() { relativeAccuracy = defaultRelativeAccuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#setFunctionValueAccuracy(double) + /** + * Set the function value accuracy. + * + * @param accuracy the accuracy. + * @throws MathException if the accuracy can't be achieved by the solver or + * is otherwise deemed unreasonable. */ public void setFunctionValueAccuracy(double accuracy) throws MathException { functionValueAccuracy = accuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#getFunctionValueAccuracy() + /** + * Get the actual function value accuracy. + * @return the accuracy */ public double getFunctionValueAccuracy() { return functionValueAccuracy; } - /* (non-Javadoc) - * @see org.apache.commons.math.UnivariateRealSolver#resetFunctionValueAccuracy() + /** + * Reset the actual function accuracy to the default. */ public void resetFunctionValueAccuracy() { functionValueAccuracy = defaultFunctionValueAccuracy; } - } diff --git a/src/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java b/src/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java index 9003ef393..12e45c1fb 100644 --- a/src/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java +++ b/src/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java @@ -53,8 +53,10 @@ */ package org.apache.commons.math.stat.distribution; +import org.apache.commons.math.MathException; import org.apache.commons.math.analysis.RootFinding; -import org.apache.commons.math.analysis.UnivariateFunction; +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.analysis.UnivariateRealSolverFactory; /** * Base class for various continuous distributions. It provides default @@ -101,22 +103,37 @@ public abstract class AbstractContinuousDistribution // by default, do simple root finding using bracketing and bisection. // subclasses can overide if there is a better method. - UnivariateFunction rootFindingFunction = new UnivariateFunction() { - public double evaluate(double x) { + UnivariateRealFunction rootFindingFunction = + new UnivariateRealFunction() { + + public double value(double x) throws MathException { return cummulativeProbability(x) - p; } + + public double firstDerivative(double x) throws MathException { + return 0; + } + + public double secondDerivative(double x) throws MathException { + return 0; + } }; - // bracket root - double[] bracket = RootFinding.bracket(rootFindingFunction, - getInitialDomain(p), getDomainLowerBound(p), - getDomainUpperBound(p)); + try { + // bracket root + double[] bracket = RootFinding.bracket(rootFindingFunction, + getInitialDomain(p), getDomainLowerBound(p), + getDomainUpperBound(p)); + + // find root + double root = UnivariateRealSolverFactory.solve( + rootFindingFunction, bracket[0], bracket[1]); - // find root - double root = RootFinding.bisection(rootFindingFunction, bracket[0], - bracket[1]); - - return root; + return root; + } catch (MathException ex) { + // this should never happen. + return Double.NaN; + } } /** diff --git a/src/test/org/apache/commons/math/analysis/BisectionSolverTest.java b/src/test/org/apache/commons/math/analysis/BisectionSolverTest.java new file mode 100644 index 000000000..0bacca012 --- /dev/null +++ b/src/test/org/apache/commons/math/analysis/BisectionSolverTest.java @@ -0,0 +1,127 @@ +/* ==================================================================== + * The Apache Software License, Version 1.1 + * + * Copyright (c) 2003 The Apache Software Foundation. All rights + * reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * 3. The end-user documentation included with the redistribution, if + * any, must include the following acknowlegement: + * "This product includes software developed by the + * Apache Software Foundation (http://www.apache.org/)." + * Alternately, this acknowlegement may appear in the software itself, + * if and wherever such third-party acknowlegements normally appear. + * + * 4. The names "The Jakarta Project", "Commons", and "Apache Software + * Foundation" must not be used to endorse or promote products derived + * from this software without prior written permission. For written + * permission, please contact apache@apache.org. + * + * 5. Products derived from this software may not be called "Apache" + * nor may "Apache" appear in their names without prior written + * permission of the Apache Software Foundation. + * + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF + * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT + * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * ==================================================================== + * + * This software consists of voluntary contributions made by many + * individuals on behalf of the Apache Software Foundation. For more + * information on the Apache Software Foundation, please see + * . + */ +package org.apache.commons.math.analysis; + +import org.apache.commons.math.MathException; + +import junit.framework.TestCase; + +/** + * + */ +public final class BisectionSolverTest extends TestCase { + /** + * + */ + public BisectionSolverTest(String name) { + super(name); + } + + /** + * + */ + public void testSinZero() throws MathException { + UnivariateRealFunction f = new SinFunction(); + double result; + + UnivariateRealSolver solver = new BisectionSolver(f); + result = solver.solve(3, 4); + assertEquals(result, Math.PI, solver.getAbsoluteAccuracy()); + + result = solver.solve(1, 4); + assertEquals(result, Math.PI, solver.getAbsoluteAccuracy()); + } + + /** + * + */ + public void testQuinticZero() throws MathException { + UnivariateRealFunction f = new QuinticFunction(); + double result; + + UnivariateRealSolver solver = new BisectionSolver(f); + result = solver.solve(-0.2, 0.2); + assertEquals(result, 0, solver.getAbsoluteAccuracy()); + + result = solver.solve(-0.1, 0.3); + assertEquals(result, 0, solver.getAbsoluteAccuracy()); + + result = solver.solve(-0.3, 0.45); + assertEquals(result, 0, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.3, 0.7); + assertEquals(result, 0.5, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.2, 0.6); + assertEquals(result, 0.5, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.05, 0.95); + assertEquals(result, 0.5, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.85, 1.25); + assertEquals(result, 1.0, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.8, 1.2); + assertEquals(result, 1.0, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.85, 1.75); + assertEquals(result, 1.0, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.55, 1.45); + assertEquals(result, 1.0, solver.getAbsoluteAccuracy()); + + result = solver.solve(0.85, 5); + assertEquals(result, 1.0, solver.getAbsoluteAccuracy()); + } +}