Mark R. Diggory 2003-07-05 18:05:31 +00:00
parent 561a4c85e7
commit 2b5b0e7a14
8 changed files with 501 additions and 172 deletions

View File

@ -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
* <http://www.apache.org/>.
*/
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;
}
}

View File

@ -64,16 +64,36 @@ import org.apache.commons.math.MathException;
* @author pietsch at apache.org * @author pietsch at apache.org
*/ */
public class BrentSolver extends UnivariateRealSolverImpl { public class BrentSolver extends UnivariateRealSolverImpl {
/**
private UnivariateRealFunction f; * Construct a solver for the given function.
* @param f function to solve.
*/
public BrentSolver(UnivariateRealFunction f) { public BrentSolver(UnivariateRealFunction f) {
super(100, 1E-6); super(f, 100, 1E-6);
this.f = f;
} }
/* (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 { public double solve(double min, double max) throws MathException {
clearResult(); clearResult();
@ -116,21 +136,8 @@ public class BrentSolver extends UnivariateRealSolverImpl {
setResult(x1, i); setResult(x1, i);
return result; 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 if (Math.abs(oldDelta) < tolerance
|| Math.abs(y0) <= Math.abs(y1)) { || Math.abs(y0) <= Math.abs(y1)) {
// System.out.println("bisection");
// Force bisection. // Force bisection.
delta = 0.5 * dx; delta = 0.5 * dx;
oldDelta = delta; oldDelta = delta;
@ -140,12 +147,10 @@ public class BrentSolver extends UnivariateRealSolverImpl {
double p1; double p1;
if (x0 == x2) { if (x0 == x2) {
// Linear interpolation. // Linear interpolation.
// System.out.println("linear");
p = dx * r3; p = dx * r3;
p1 = 1.0 - r3; p1 = 1.0 - r3;
} else { } else {
// Inverse quadratic interpolation. // Inverse quadratic interpolation.
// System.out.println("invers quad");
double r1 = y0 / y2; double r1 = y0 / y2;
double r2 = y1 / y2; double r2 = y1 / y2;
p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0)); p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
@ -156,14 +161,11 @@ public class BrentSolver extends UnivariateRealSolverImpl {
} else { } else {
p = -p; 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) if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1)
|| p >= Math.abs(0.5 * oldDelta * p1)) { || p >= Math.abs(0.5 * oldDelta * p1)) {
// Inverse quadratic interpolation gives a value // Inverse quadratic interpolation gives a value
// in the wrong direction, or progress is slow. // in the wrong direction, or progress is slow.
// Fall back to bisection. // Fall back to bisection.
// System.out.println("bisection fallback");
delta = 0.5 * dx; delta = 0.5 * dx;
oldDelta = delta; oldDelta = delta;
} else { } else {
@ -178,9 +180,9 @@ public class BrentSolver extends UnivariateRealSolverImpl {
if (Math.abs(delta) > tolerance) { if (Math.abs(delta) > tolerance) {
x1 = x1 + delta; x1 = x1 + delta;
} else if (dx > 0.0) { } else if (dx > 0.0) {
x1 = x1 + 0.5*tolerance; x1 = x1 + 0.5 * tolerance;
} else if (dx <= 0.0) { } else if (dx <= 0.0) {
x1 = x1 - 0.5*tolerance; x1 = x1 - 0.5 * tolerance;
} }
y1 = f.value(x1); y1 = f.value(x1);
if ((y1 > 0) == (y2 > 0)) { if ((y1 > 0) == (y2 > 0)) {
@ -193,5 +195,4 @@ public class BrentSolver extends UnivariateRealSolverImpl {
} }
throw new MathException("Maximal iteration number exceeded."); throw new MathException("Maximal iteration number exceeded.");
} }
} }

View File

@ -53,6 +53,8 @@
*/ */
package org.apache.commons.math.analysis; package org.apache.commons.math.analysis;
import org.apache.commons.math.MathException;
/** /**
* Utility class comprised of root finding techniques. * Utility class comprised of root finding techniques.
* *
@ -77,14 +79,17 @@ public class RootFinding {
* @param function the function * @param function the function
* @param initial midpoint of the returned range. * @param initial midpoint of the returned range.
* @param lowerBound for numerical safety, a never is less than this value. * @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}. * @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 initial,
double lowerBound, double lowerBound,
double upperBound) { double upperBound) throws MathException {
return bracket( function, initial, lowerBound, upperBound, Integer.MAX_VALUE ) ; return bracket( function, initial, lowerBound, upperBound,
Integer.MAX_VALUE ) ;
} }
/** /**
@ -95,15 +100,18 @@ public class RootFinding {
* @param function the function * @param function the function
* @param initial midpoint of the returned range. * @param initial midpoint of the returned range.
* @param lowerBound for numerical safety, a never is less than this value. * @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
* @param maximumIterations to guard against infinite looping, maximum number of iterations to perform * value.
* @param maximumIterations to guard against infinite looping, maximum
* number of iterations to perform
* @return a two element array holding {a, b}. * @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 initial,
double lowerBound, double lowerBound,
double upperBound, double upperBound,
int maximumIterations) { int maximumIterations) throws MathException {
double a = initial; double a = initial;
double b = initial; double b = initial;
double fa; double fa;
@ -113,52 +121,11 @@ public class RootFinding {
do { do {
a = Math.max(a - 1.0, lowerBound); a = Math.max(a - 1.0, lowerBound);
b = Math.min(b + 1.0, upperBound); b = Math.min(b + 1.0, upperBound);
fa = function.evaluate(a); fa = function.value(a);
fb = function.evaluate(b); fb = function.value(b);
numIterations += 1 ; numIterations += 1 ;
} while ( (fa * fb > 0.0) && ( numIterations < maximumIterations ) ); } while ( (fa * fb > 0.0) && ( numIterations < maximumIterations ) );
return new double[]{a, b}; 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;
}
} }

View File

@ -66,16 +66,36 @@ import org.apache.commons.math.MathException;
* @author pietsch at apache.org * @author pietsch at apache.org
*/ */
public class SecantSolver extends UnivariateRealSolverImpl { public class SecantSolver extends UnivariateRealSolverImpl {
/**
private UnivariateRealFunction f; * Construct a solver for the given function.
* @param f function to solve.
*/
public SecantSolver(UnivariateRealFunction f) { public SecantSolver(UnivariateRealFunction f) {
super(100, 1E-6); super(f, 100, 1E-6);
this.f = f;
} }
/* (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 { public double solve(double min, double max) throws MathException {
clearResult(); clearResult();
@ -88,7 +108,7 @@ public class SecantSolver extends UnivariateRealSolverImpl {
double x1 = max; double x1 = max;
double y0 = f.value(x0); double y0 = f.value(x0);
double y1 = f.value(x1); double y1 = f.value(x1);
if ((y0>0)== (y1>0)) { if ((y0 > 0) == (y1 > 0)) {
throw new MathException("Interval doesn't bracket a zero."); throw new MathException("Interval doesn't bracket a zero.");
} }
double x2 = x0; double x2 = x0;
@ -117,21 +137,19 @@ public class SecantSolver extends UnivariateRealSolverImpl {
if (Math.abs(y1) > Math.abs(y0)) { if (Math.abs(y1) > Math.abs(y0)) {
// Function value increased in last iteration. Force bisection. // Function value increased in last iteration. Force bisection.
delta = 0.5 * oldDelta; delta = 0.5 * oldDelta;
// System.out.println("Forced Bisection");
} else { } else {
delta = (x0 - x1) / (1 - y0 / y1); delta = (x0 - x1) / (1 - y0 / y1);
// System.out.println("delta=" + delta + " olddelta=" + oldDelta);
if (delta / oldDelta > 1) { 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; delta = 0.5 * oldDelta;
// System.out.println("Fallback Bisection");
} }
} }
x0 = x1; x0 = x1;
y0 = y1; y0 = y1;
x1 = x1 + delta; x1 = x1 + delta;
y1 = f.value(x1); y1 = f.value(x1);
if ((y1>0) == (y2>0)) { if ((y1 > 0) == (y2 > 0)) {
// New bracket is (x0,x1). // New bracket is (x0,x1).
x2 = x0; x2 = x0;
y2 = y0; y2 = y0;

View File

@ -66,23 +66,27 @@ public interface UnivariateRealSolver {
/** /**
* Set the upper limit for the number of iterations. * Set the upper limit for the number of iterations.
*
* Usually a high iteration count indicates convergence problems. However, * Usually a high iteration count indicates convergence problems. However,
* the "reasonable value" varies widely for different solvers, users are * 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. * An exception will be thrown if the number is exceeded.
* *
* @param count * @param count maximum number of iterations
*/ */
public void setMaximalIterationCount(int count); public void setMaximalIterationCount(int count);
/** /**
* Get the upper limit for the number of iterations. * Get the upper limit for the number of iterations.
*
* @return the actual upper limit * @return the actual upper limit
*/ */
public int getMaximalIterationCount(); public int getMaximalIterationCount();
/** /**
* Reset the upper limit for the number of iterations to the default. * Reset the upper limit for the number of iterations to the default.
*
* The default value is supplied by the solver implementation. * The default value is supplied by the solver implementation.
* *
* @see #setMaximalIterationCount(int) * @see #setMaximalIterationCount(int)
@ -91,41 +95,50 @@ public interface UnivariateRealSolver {
/** /**
* Set the absolute accuracy. * Set the absolute accuracy.
*
* The default is usually choosen so taht roots in the interval * 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 * -10..-0.1 and +0.1..+10 can be found wit a reasonable accuracy. If the
* absolute value of your roots is of much smaller magnitude, set this to a smaller * expected absolute value of your roots is of much smaller magnitude, set
* value. * this to a smaller value.
* Solvers are advised to do a plausibility check with the relative accuracy, but *
* clients should not rely on this. * Solvers are advised to do a plausibility check with the relative
* accuracy, but clients should not rely on this.
*
* @param accuracy the accuracy. * @param accuracy the accuracy.
* @throws MathException if the accuracy can't be achieved by the solver or is * @throws MathException if the accuracy can't be achieved by the solver or
* otherwise deemed unreasonable. * is otherwise deemed unreasonable.
*/ */
public void setAbsoluteAccuracy(double accuracy) throws MathException; public void setAbsoluteAccuracy(double accuracy) throws MathException;
/** /**
* Get the actual absolute accuracy. * Get the actual absolute accuracy.
*
* @return the accuracy * @return the accuracy
*/ */
public double getAbsoluteAccuracy(); public double getAbsoluteAccuracy();
/** /**
* Reset the absolute accuracy to the default. * Reset the absolute accuracy to the default.
*
* The default value is provided by the solver implementation. * The default value is provided by the solver implementation.
*/ */
public void resetAbsoluteAccuracy(); public void resetAbsoluteAccuracy();
/** /**
* Set the relative accuracy. * 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. * This is used to stop iterations if the absolute accuracy can't be
* If this should be the primary criterium for convergence rather then a safety * achieved due to large values or short mantissa length.
* measure, set the absolute accuracy to a ridiculously small value, like 1E-1000. *
* 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. * @param accuracy the relative accuracy.
* @throws MathException if the accuracy can't be achieved by the solver or is * @throws MathException if the accuracy can't be achieved by the solver or
* otherwise deemed unreasonable. * is otherwise deemed unreasonable.
*/ */
public void setRelativeAccuracy(double Accuracy) throws MathException; public void setRelativeAccuracy(double accuracy) throws MathException;
/** /**
* Get the actual relative accuracy. * Get the actual relative accuracy.
@ -141,14 +154,18 @@ public interface UnivariateRealSolver {
/** /**
* Set the function value accuracy. * Set the function value accuracy.
*
* This is used to determine whan an evaluated function value or some other * This is used to determine whan an evaluated function value or some other
* value which is used as divisor is zero. * 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. * @param accuracy the accuracy.
* @throws MathException if the accuracy can't be achieved by the solver or is * @throws MathException if the accuracy can't be achieved by the solver or
* otherwise deemed unreasonable. * is otherwise deemed unreasonable.
*/ */
public void setFunctionValueAccuracy(double Accuracy) throws MathException; public void setFunctionValueAccuracy(double accuracy) throws MathException;
/** /**
* Get the actual function value accuracy. * Get the actual function value accuracy.

View File

@ -65,24 +65,57 @@ import org.apache.commons.math.MathException;
public abstract class UnivariateRealSolverImpl public abstract class UnivariateRealSolverImpl
implements UnivariateRealSolver { implements UnivariateRealSolver {
/** Maximum absolute error. */
protected double absoluteAccuracy; protected double absoluteAccuracy;
/** Maximum relative error. */
protected double relativeAccuracy; protected double relativeAccuracy;
/** Maximum error of function. */
protected double functionValueAccuracy; protected double functionValueAccuracy;
/** Maximum number of iterations. */
protected int maximalIterationCount; protected int maximalIterationCount;
/** Default maximum absolute error. */
protected double defaultAbsoluteAccuracy; protected double defaultAbsoluteAccuracy;
/** Default maximum relative error. */
protected double defaultRelativeAccuracy; protected double defaultRelativeAccuracy;
/** Default maximum error of function. */
protected double defaultFunctionValueAccuracy; protected double defaultFunctionValueAccuracy;
/** Default maximum number of iterations. */
protected int defaultMaximalIterationCount; protected int defaultMaximalIterationCount;
/** Indicates where a root has been computed. */
protected boolean resultComputed = false; protected boolean resultComputed = false;
/** The last computed root. */
protected double result; protected double result;
// Mainly for test framework. // Mainly for test framework.
/** The last iteration count. */
protected int iterationCount; 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( protected UnivariateRealSolverImpl(
UnivariateRealFunction f,
int defaultMaximalIterationCount, int defaultMaximalIterationCount,
double defaultAbsoluteAccuracy) { double defaultAbsoluteAccuracy) {
super();
this.f = f;
this.defaultAbsoluteAccuracy = defaultAbsoluteAccuracy; this.defaultAbsoluteAccuracy = defaultAbsoluteAccuracy;
this.defaultRelativeAccuracy = 1E-14; this.defaultRelativeAccuracy = 1E-14;
this.defaultFunctionValueAccuracy = 1E-15; this.defaultFunctionValueAccuracy = 1E-15;
@ -93,48 +126,39 @@ public abstract class UnivariateRealSolverImpl
this.maximalIterationCount = defaultMaximalIterationCount; this.maximalIterationCount = defaultMaximalIterationCount;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#solve(double, double) * Access the last computed root.
*/ * @return the last computed root.
public double solve(double min, double max) throws MathException { * @throws MathException if no root has been computed.
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()
*/ */
public double getResult() throws MathException { public double getResult() throws MathException {
if (resultComputed) { if (resultComputed) {
return result; return result;
} else { } else {
// TODO: could this be an IllegalStateException instead?
throw new MathException("No result available"); 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 { public int getIterationCount() throws MathException {
if (resultComputed) { if (resultComputed) {
return iterationCount; return iterationCount;
} else { } else {
// TODO: could this be an IllegalStateException instead?
throw new MathException("No result available"); throw new MathException("No result available");
} }
} }
/* /**
* Convenience function for implementations. * Convenience function for implementations.
* @param result the result to set * @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) { protected final void setResult(double result, int iterationCount) {
this.result = result; this.result = result;
@ -142,97 +166,116 @@ public abstract class UnivariateRealSolverImpl
this.resultComputed = true; this.resultComputed = true;
} }
/* /**
* Convenience function for implementations. * Convenience function for implementations.
*/ */
protected final void clearResult() { protected final void clearResult() {
this.resultComputed = false; 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) public void setAbsoluteAccuracy(double accuracy)
throws MathException { throws MathException {
absoluteAccuracy = accuracy; absoluteAccuracy = accuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#getAccuracy() * Get the actual absolute accuracy.
*
* @return the accuracy
*/ */
public double getAbsoluteAccuracy() { public double getAbsoluteAccuracy() {
return absoluteAccuracy; return absoluteAccuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#resetAbsoluteAccuracy() * Reset the absolute accuracy to the default.
*/ */
public void resetAbsoluteAccuracy() { public void resetAbsoluteAccuracy() {
absoluteAccuracy = defaultAbsoluteAccuracy; 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) { public void setMaximalIterationCount(int count) {
maximalIterationCount = 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() { public int getMaximalIterationCount() {
return maximalIterationCount; 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() { public void resetMaximalIterationCount() {
maximalIterationCount = defaultMaximalIterationCount; 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 { public void setRelativeAccuracy(double accuracy) throws MathException {
relativeAccuracy = accuracy; relativeAccuracy = accuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#getRelativeAccuracy() * Get the actual relative accuracy.
* @return the accuracy
*/ */
public double getRelativeAccuracy() { public double getRelativeAccuracy() {
return relativeAccuracy; return relativeAccuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#resetRelativeAccuracy() * Reset the relative accuracy to the default.
*/ */
public void resetRelativeAccuracy() { public void resetRelativeAccuracy() {
relativeAccuracy = defaultRelativeAccuracy; 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) public void setFunctionValueAccuracy(double accuracy)
throws MathException { throws MathException {
functionValueAccuracy = accuracy; functionValueAccuracy = accuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#getFunctionValueAccuracy() * Get the actual function value accuracy.
* @return the accuracy
*/ */
public double getFunctionValueAccuracy() { public double getFunctionValueAccuracy() {
return functionValueAccuracy; return functionValueAccuracy;
} }
/* (non-Javadoc) /**
* @see org.apache.commons.math.UnivariateRealSolver#resetFunctionValueAccuracy() * Reset the actual function accuracy to the default.
*/ */
public void resetFunctionValueAccuracy() { public void resetFunctionValueAccuracy() {
functionValueAccuracy = defaultFunctionValueAccuracy; functionValueAccuracy = defaultFunctionValueAccuracy;
} }
} }

View File

@ -53,8 +53,10 @@
*/ */
package org.apache.commons.math.stat.distribution; 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.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 * 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. // by default, do simple root finding using bracketing and bisection.
// subclasses can overide if there is a better method. // subclasses can overide if there is a better method.
UnivariateFunction rootFindingFunction = new UnivariateFunction() { UnivariateRealFunction rootFindingFunction =
public double evaluate(double x) { new UnivariateRealFunction() {
public double value(double x) throws MathException {
return cummulativeProbability(x) - p; return cummulativeProbability(x) - p;
} }
public double firstDerivative(double x) throws MathException {
return 0;
}
public double secondDerivative(double x) throws MathException {
return 0;
}
}; };
try {
// bracket root // bracket root
double[] bracket = RootFinding.bracket(rootFindingFunction, double[] bracket = RootFinding.bracket(rootFindingFunction,
getInitialDomain(p), getDomainLowerBound(p), getInitialDomain(p), getDomainLowerBound(p),
getDomainUpperBound(p)); getDomainUpperBound(p));
// find root // find root
double root = RootFinding.bisection(rootFindingFunction, bracket[0], double root = UnivariateRealSolverFactory.solve(
bracket[1]); rootFindingFunction, bracket[0], bracket[1]);
return root; return root;
} catch (MathException ex) {
// this should never happen.
return Double.NaN;
}
} }
/** /**

View File

@ -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
* <http://www.apache.org/>.
*/
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());
}
}