Added numerical integration classes contributed by Xiaogang Zhang.

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@232588 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz 2005-08-14 06:52:31 +00:00
parent bd34147ab1
commit 11796d4ed1
9 changed files with 1089 additions and 1 deletions

View File

@ -0,0 +1,108 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
/**
* Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
* Romberg Algorithm</a> for integrating of real univariate functions. For
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
* chapter 3.
* <p>
* Romberg integration employs k successvie refinements of the trapezoid
* rule to remove error terms less than order O(N^(-2k)). Simpson's rule
* is a special case of k = 2.
*
* @version $Revision$ $Date$
*/
public class RombergIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = -1058849527738180243L;
/**
* Construct an integrator for the given function.
*
* @param f function to solve
*/
public RombergIntegrator(UnivariateRealFunction f) {
super(f, 32);
}
/**
* Integrate the function in the given interval.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1, j, m = maximalIterationCount + 1;
// Array strcture here can be improved for better space
// efficiency because only the lower triangle is used.
double r, t[][] = new double[m][m], s, olds;
clearResult();
verifyInterval(min, max);
verifyIterationCount();
TrapezoidIntegrator qtrap = new TrapezoidIntegrator(this.f);
t[0][0] = qtrap.stage(min, max, 0);
olds = t[0][0];
while (i <= maximalIterationCount) {
t[i][0] = qtrap.stage(min, max, i);
for (j = 1; j <= i; j++) {
// Richardson extrapolation coefficient
r = (1L << (2 * j)) -1;
t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
}
s = t[i][i];
if (i >= minimalIterationCount) {
if (Math.abs(s - olds) <= Math.abs(relativeAccuracy * olds)) {
setResult(s, i);
return result;
}
}
olds = s;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
/**
* Verifies that the iteration limits are valid and within the range.
*
* @throws IllegalArgumentException if not
*/
protected void verifyIterationCount() throws IllegalArgumentException {
super.verifyIterationCount();
// at most 32 bisection refinements due to higher order divider
if (maximalIterationCount > 32) {
throw new IllegalArgumentException
("Iteration upper limit out of [0, 32] range: " +
maximalIterationCount);
}
}
}

View File

@ -0,0 +1,107 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
/**
* Implements the <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
* Simpson's Rule</a> for integrating of real univariate functions. For
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
* chapter 3.
* <p>
* This implementation employs basic trapezoid rule as building blocks to
* calculate the Simpson's rule of alternating 2/3 and 4/3.
*
* @version $Revision$ $Date$
*/
public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = 3405465123320678216L;
/**
* Construct an integrator for the given function.
*
* @param f function to solve
*/
public SimpsonIntegrator(UnivariateRealFunction f) {
super(f, 64);
}
/**
* Integrate the function in the given interval.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1;
double s, olds, t, oldt;
clearResult();
verifyInterval(min, max);
verifyIterationCount();
TrapezoidIntegrator qtrap = new TrapezoidIntegrator(this.f);
if (minimalIterationCount == 1) {
s = (4 * qtrap.stage(min, max, 1) - qtrap.stage(min, max, 0)) / 3.0;
setResult(s, 1);
return result;
}
// Simpson's rule requires at least two trapezoid stages.
olds = 0;
oldt = qtrap.stage(min, max, 0);
while (i <= maximalIterationCount) {
t = qtrap.stage(min, max, i);
s = (4 * t - oldt) / 3.0;
if (i >= minimalIterationCount) {
if (Math.abs(s - olds) <= Math.abs(relativeAccuracy * olds)) {
setResult(s, i);
return result;
}
}
olds = s;
oldt = t;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
/**
* Verifies that the iteration limits are valid and within the range.
*
* @throws IllegalArgumentException if not
*/
protected void verifyIterationCount() throws IllegalArgumentException {
super.verifyIterationCount();
// at most 64 bisection refinements
if (maximalIterationCount > 64) {
throw new IllegalArgumentException
("Iteration upper limit out of [0, 64] range: " +
maximalIterationCount);
}
}
}

View File

@ -0,0 +1,138 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
/**
* Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html">
* Trapezoidal Rule</a> for integrating of real univariate functions. For
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
* chapter 3.
* <p>
* The function should be integrable.
*
* @version $Revision$ $Date$
*/
public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = 4978222553983172543L;
/** intermediate result */
private double s;
/**
* Construct an integrator for the given function.
*
* @param f function to solve
*/
public TrapezoidIntegrator(UnivariateRealFunction f) {
super(f, 64);
}
/**
* Compute the n-th stage integral of trapezoid rule. This function
* should only be called by API <code>integrate()</code> in the package.
* To save time it does not verify arguments - caller does.
* <p>
* The interval is divided equally into 2^n sections rather than an
* arbitrary m sections because this configuration can best utilize the
* alrealy computed values.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param n the stage of 1/2 refinement, n = 0 is no refinement
* @return the value of n-th stage integral
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
*/
double stage(double min, double max, int n) throws
FunctionEvaluationException {
long i, np;
double x, spacing, sum = 0;
if (n == 0) {
s = 0.5 * (max - min) * (f.value(min) + f.value(max));
return s;
} else {
np = 1L << (n-1); // number of new points in this stage
spacing = (max - min) / np; // spacing between adjacent new points
x = min + 0.5 * spacing; // the first new point
for (i = 0; i < np; i++) {
sum += f.value(x);
x += spacing;
}
// add the new sum to previously calculated result
s = 0.5 * (s + sum * spacing);
return s;
}
}
/**
* Integrate the function in the given interval.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1;
double t, oldt;
clearResult();
verifyInterval(min, max);
verifyIterationCount();
oldt = stage(min, max, 0);
while (i <= maximalIterationCount) {
t = stage(min, max, i);
if (i >= minimalIterationCount) {
if (Math.abs(t - oldt) <= Math.abs(relativeAccuracy * oldt)) {
setResult(t, i);
return result;
}
}
oldt = t;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
/**
* Verifies that the iteration limits are valid and within the range.
*
* @throws IllegalArgumentException if not
*/
protected void verifyIterationCount() throws IllegalArgumentException {
super.verifyIterationCount();
// at most 64 bisection refinements
if (maximalIterationCount > 64) {
throw new IllegalArgumentException
("Iteration upper limit out of [0, 64] range: " +
maximalIterationCount);
}
}
}

View File

@ -0,0 +1,154 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
/**
* Interface for univariate real integration algorithms.
*
* @version $Revision$ $Date$
*/
public interface UnivariateRealIntegrator {
/**
* Set the upper limit for the number of iterations.
* <p>
* Usually a high iteration count indicates convergence problem. However,
* the "reasonable value" varies widely for different cases. Users are
* advised to use the default value.
* <p>
* A <code>ConvergenceException</code> will be thrown if this number
* is exceeded.
*
* @param count maximum number of iterations
*/
void setMaximalIterationCount(int count);
/**
* Get the upper limit for the number of iterations.
*
* @return the actual upper limit
*/
int getMaximalIterationCount();
/**
* Reset the upper limit for the number of iterations to the default.
* <p>
* The default value is supplied by the implementation.
*
* @see #setMaximalIterationCount(int)
*/
void resetMaximalIterationCount();
/**
* Set the lower limit for the number of iterations.
* <p>
* Minimal iteration is needed to avoid false early convergence, e.g.
* the sample points happen to be zeroes of the function. Users can
* use the default value or choose one that they see as appropriate.
* <p>
* A <code>ConvergenceException</code> will be thrown if this number
* is not met.
*
* @param count minimum number of iterations
*/
void setMinimalIterationCount(int count);
/**
* Get the lower limit for the number of iterations.
*
* @return the actual lower limit
*/
int getMinimalIterationCount();
/**
* Reset the lower limit for the number of iterations to the default.
* <p>
* The default value is supplied by the implementation.
*
* @see #setMinimalIterationCount(int)
*/
void resetMinimalIterationCount();
/**
* Set the relative accuracy.
* <p>
* This is used to stop iterations.
*
* @param accuracy the relative accuracy
* @throws IllegalArgumentException if the accuracy can't be achieved
* or is otherwise deemed unreasonable
*/
void setRelativeAccuracy(double accuracy);
/**
* Get the actual relative accuracy.
*
* @return the accuracy
*/
double getRelativeAccuracy();
/**
* Reset the relative accuracy to the default.
* <p>
* The default value is provided by the implementation.
*
* @see #setRelativeAccuracy(double)
*/
void resetRelativeAccuracy();
/**
* Integrate the function in the given interval.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min > max or the endpoints do not
* satisfy the requirements specified by the integrator
*/
double integrate(double min, double max) throws ConvergenceException,
FunctionEvaluationException, IllegalArgumentException;
/**
* Get the result of the last run of the integrator.
*
* @return the last result
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed
*/
double getResult() throws IllegalStateException;
/**
* Get the number of iterations in the last run of the integrator.
* <p>
* This is mainly meant for testing purposes. It may occasionally
* help track down performance problems: if the iteration count
* is notoriously high, check whether the function is evaluated
* properly, and whether another integrator is more amenable to the
* problem.
*
* @return the last iteration count
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed
*/
int getIterationCount() throws IllegalStateException;
}

View File

@ -0,0 +1,254 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import java.io.Serializable;
/**
* Provide a default implementation for several generic functions.
*
* @version $Revision$ $Date$
*/
public abstract class UnivariateRealIntegratorImpl implements
UnivariateRealIntegrator, Serializable {
/** serializable version identifier */
static final long serialVersionUID = -3365294665201465048L;
/** maximum relative error */
protected double relativeAccuracy;
/** maximum number of iterations */
protected int maximalIterationCount;
/** minimum number of iterations */
protected int minimalIterationCount;
/** default maximum relative error */
protected double defaultRelativeAccuracy;
/** default maximum number of iterations */
protected int defaultMaximalIterationCount;
/** default minimum number of iterations */
protected int defaultMinimalIterationCount;
/** indicates whether an integral has been computed */
protected boolean resultComputed = false;
/** the last computed integral */
protected double result;
/** the last iteration count */
protected int iterationCount;
/** the integrand function */
protected UnivariateRealFunction f;
/**
* Construct an integrator with given iteration count and accuracy.
*
* @param f the integrand function
* @param defaultMaximalIterationCount maximum number of iterations
* @throws IllegalArgumentException if f is null or the iteration
* limits are not valid
*/
protected UnivariateRealIntegratorImpl(
UnivariateRealFunction f,
int defaultMaximalIterationCount) throws IllegalArgumentException {
if (f == null) {
throw new IllegalArgumentException("Function can not be null.");
}
this.f = f;
// parameters that may depend on algorithm
this.defaultMaximalIterationCount = defaultMaximalIterationCount;
this.maximalIterationCount = defaultMaximalIterationCount;
// parameters that are problem specific
this.defaultRelativeAccuracy = 1E-6;
this.relativeAccuracy = defaultRelativeAccuracy;
this.defaultMinimalIterationCount = 3;
this.minimalIterationCount = defaultMinimalIterationCount;
verifyIterationCount();
}
/**
* Access the last computed integral.
*
* @return the last computed integral
* @throws IllegalStateException if no integral has been computed
*/
public double getResult() throws IllegalStateException {
if (resultComputed) {
return result;
} else {
throw new IllegalStateException("No result available.");
}
}
/**
* Access the last iteration count.
*
* @return the last iteration count
* @throws IllegalStateException if no integral has been computed
*/
public int getIterationCount() throws IllegalStateException {
if (resultComputed) {
return iterationCount;
} else {
throw new IllegalStateException("No result available.");
}
}
/**
* Convenience function for implementations.
*
* @param result the result to set
* @param iterationCount the iteration count to set
*/
protected final void setResult(double result, int iterationCount) {
this.result = result;
this.iterationCount = iterationCount;
this.resultComputed = true;
}
/**
* Convenience function for implementations.
*/
protected final void clearResult() {
this.resultComputed = false;
}
/**
* Set the upper limit for the number of iterations.
*
* @param count maximum number of iterations
*/
public void setMaximalIterationCount(int count) {
maximalIterationCount = count;
}
/**
* Get the upper limit for the number of iterations.
*
* @return the actual upper limit
*/
public int getMaximalIterationCount() {
return maximalIterationCount;
}
/**
* Reset the upper limit for the number of iterations to the default.
*/
public void resetMaximalIterationCount() {
maximalIterationCount = defaultMaximalIterationCount;
}
/**
* Set the lower limit for the number of iterations.
*
* @param count minimum number of iterations
*/
public void setMinimalIterationCount(int count) {
minimalIterationCount = count;
}
/**
* Get the lower limit for the number of iterations.
*
* @return the actual lower limit
*/
public int getMinimalIterationCount() {
return minimalIterationCount;
}
/**
* Reset the lower limit for the number of iterations to the default.
*/
public void resetMinimalIterationCount() {
minimalIterationCount = defaultMinimalIterationCount;
}
/**
* Set the relative accuracy.
*
* @param accuracy the relative accuracy
* @throws IllegalArgumentException if the accuracy can't be achieved by
* the integrator or is otherwise deemed unreasonable
*/
public void setRelativeAccuracy(double accuracy) {
relativeAccuracy = accuracy;
}
/**
* Get the actual relative accuracy.
*
* @return the accuracy
*/
public double getRelativeAccuracy() {
return relativeAccuracy;
}
/**
* Reset the relative accuracy to the default.
*/
public void resetRelativeAccuracy() {
relativeAccuracy = defaultRelativeAccuracy;
}
/**
* Returns true if the arguments form a (strictly) increasing sequence
*
* @param start first number
* @param mid second number
* @param end third number
* @return true if the arguments form an increasing sequence
*/
protected boolean isSequence(double start, double mid, double end) {
return (start < mid) && (mid < end);
}
/**
* Verifies that the endpoints specify an interval.
*
* @param lower lower endpoint
* @param upper upper endpoint
* @throws IllegalArgumentException if not interval
*/
protected void verifyInterval(double lower, double upper) throws
IllegalArgumentException {
if (lower >= upper) {
throw new IllegalArgumentException
("Endpoints do not specify an interval: [" + lower +
", " + upper + "]");
}
}
/**
* Verifies that the upper and lower limits of iterations are valid.
*
* @throws IllegalArgumentException if not valid
*/
protected void verifyIterationCount() throws IllegalArgumentException {
if (!isSequence(0, minimalIterationCount, maximalIterationCount+1)) {
throw new IllegalArgumentException
("Invalid iteration limits: min=" + minimalIterationCount +
" max=" + maximalIterationCount);
}
}
}

View File

@ -0,0 +1,108 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Romberg integrator.
* <p>
* Romberg algorithm is very fast for good behavior integrand. Test runs
* show that for a default relative accuracy of 1E-6, it generally takes
* takes less than 5 iterations for the integral to converge.
*
* @version $Revision$ $Date$
*/
public final class RombergIntegratorTest extends TestCase {
/**
* Test of integrator for the sine function.
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new RombergIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = Math.PI; expected = 2;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -Math.PI/3; max = 0; expected = -0.5;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of integrator for the quintic function.
*/
public void testQuinticFunction() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealIntegrator integrator = new RombergIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = 1; expected = -1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = 0; max = 0.5; expected = 11.0/768;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the integrator.
*/
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new RombergIntegrator(f);
try {
// bad interval
integrator.integrate(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(5);
integrator.setMaximalIterationCount(4);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(10);
integrator.setMaximalIterationCount(50);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

View File

@ -0,0 +1,107 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Simpson integrator.
* <p>
* Test runs show that for a default relative accuracy of 1E-6, it
* generally takes 5 to 10 iterations for the integral to converge.
*
* @version $Revision$ $Date$
*/
public final class SimpsonIntegratorTest extends TestCase {
/**
* Test of integrator for the sine function.
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new SimpsonIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = Math.PI; expected = 2;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -Math.PI/3; max = 0; expected = -0.5;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of integrator for the quintic function.
*/
public void testQuinticFunction() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealIntegrator integrator = new SimpsonIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = 1; expected = -1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = 0; max = 0.5; expected = 11.0/768;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the integrator.
*/
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new SimpsonIntegrator(f);
try {
// bad interval
integrator.integrate(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(5);
integrator.setMaximalIterationCount(4);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(10);
integrator.setMaximalIterationCount(99);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

View File

@ -0,0 +1,107 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for trapezoid integrator.
* <p>
* Test runs show that for a default relative accuracy of 1E-6, it
* generally takes 10 to 15 iterations for the integral to converge.
*
* @version $Revision$ $Date$
*/
public final class TrapezoidIntegratorTest extends TestCase {
/**
* Test of integrator for the sine function.
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new TrapezoidIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = Math.PI; expected = 2;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -Math.PI/3; max = 0; expected = -0.5;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of integrator for the quintic function.
*/
public void testQuinticFunction() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealIntegrator integrator = new TrapezoidIntegrator(f);
double min, max, expected, result, tolerance;
min = 0; max = 1; expected = -1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = 0; max = 0.5; expected = 11.0/768;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
tolerance = Math.abs(expected * integrator.getRelativeAccuracy());
result = integrator.integrate(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the integrator.
*/
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealIntegrator integrator = new TrapezoidIntegrator(f);
try {
// bad interval
integrator.integrate(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(5);
integrator.setMaximalIterationCount(4);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad iteration limits
integrator.setMinimalIterationCount(10);
integrator.setMaximalIterationCount(99);
integrator.integrate(-1, 1);
fail("Expecting IllegalArgumentException - bad iteration limits");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

View File

@ -38,7 +38,12 @@ The <action> type attribute can be add,update,fix,remove.
Commons Math Release Notes</title>
</properties>
<body>
<release version="1.1-RC1" date="TBD"
<release version="1.2" date="TBD">
<action dev="psteitz" type="add" due-to="Xiaogang Zhang">
Added Trapeziod, Simpson, Romberg numerical integration.
</action>
</release>
<release version="1.1" date="TBD"
description="This is a maintenance release containing bug fixes and enhancements.
All API changes are binary compatible with version 1.0. The enhancements
include some new probability distributions, a Fraction class, new matrix