Added midpoint integration method.
Patch contributed by Oleksandr Kornieiev. JIRA: MATH-967 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1488914 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
d270055e87
commit
9e62190606
3
pom.xml
3
pom.xml
|
|
@ -216,6 +216,9 @@
|
|||
<contributor>
|
||||
<name>Eugene Kirpichov</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Oleksandr Kornieiev</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Piotr Kochanski</name>
|
||||
</contributor>
|
||||
|
|
|
|||
|
|
@ -51,6 +51,9 @@ If the output is not quite correct, check for invisible trailing spaces!
|
|||
</properties>
|
||||
<body>
|
||||
<release version="x.y" date="TBD" description="TBD">
|
||||
<action dev="luc" type="add" issue="MATH-967" due-to="Oleksandr Kornieiev">
|
||||
Added midpoint integration method.
|
||||
</action>
|
||||
<action dev="luc" type="fix" issue="MATH-988" due-to="Andreas Huber">
|
||||
Fixed NullPointerException in 2D and 3D sub-line intersections.
|
||||
</action>
|
||||
|
|
|
|||
|
|
@ -0,0 +1,165 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math3.analysis.integration;
|
||||
|
||||
import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
||||
import org.apache.commons.math3.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math3.exception.TooManyEvaluationsException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
|
||||
/**
|
||||
* Implements the <a href="http://en.wikipedia.org/wiki/Midpoint_method">
|
||||
* Midpoint Rule</a> for integration of real univariate functions. For
|
||||
* reference, see <b>Numerical Mathematics</b>, ISBN 0387989595,
|
||||
* chapter 9.2.
|
||||
* <p>
|
||||
* The function should be integrable.</p>
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 3.3
|
||||
*/
|
||||
public class MidPointIntegrator extends BaseAbstractUnivariateIntegrator {
|
||||
|
||||
/** Maximum number of iterations for midpoint. */
|
||||
public static final int MIDPOINT_MAX_ITERATIONS_COUNT = 64;
|
||||
|
||||
/** Intermediate result. */
|
||||
private double s;
|
||||
|
||||
/**
|
||||
* Build a midpoint integrator with given accuracies and iterations counts.
|
||||
* @param relativeAccuracy relative accuracy of the result
|
||||
* @param absoluteAccuracy absolute accuracy of the result
|
||||
* @param minimalIterationCount minimum number of iterations
|
||||
* @param maximalIterationCount maximum number of iterations
|
||||
* (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
|
||||
* @exception NotStrictlyPositiveException if minimal number of iterations
|
||||
* is not strictly positive
|
||||
* @exception NumberIsTooSmallException if maximal number of iterations
|
||||
* is lesser than or equal to the minimal number of iterations
|
||||
* @exception NumberIsTooLargeException if maximal number of iterations
|
||||
* is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
|
||||
*/
|
||||
public MidPointIntegrator(final double relativeAccuracy,
|
||||
final double absoluteAccuracy,
|
||||
final int minimalIterationCount,
|
||||
final int maximalIterationCount)
|
||||
throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
|
||||
super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
|
||||
if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
|
||||
throw new NumberIsTooLargeException(maximalIterationCount,
|
||||
MIDPOINT_MAX_ITERATIONS_COUNT, false);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Build a midpoint integrator with given iteration counts.
|
||||
* @param minimalIterationCount minimum number of iterations
|
||||
* @param maximalIterationCount maximum number of iterations
|
||||
* (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
|
||||
* @exception NotStrictlyPositiveException if minimal number of iterations
|
||||
* is not strictly positive
|
||||
* @exception NumberIsTooSmallException if maximal number of iterations
|
||||
* is lesser than or equal to the minimal number of iterations
|
||||
* @exception NumberIsTooLargeException if maximal number of iterations
|
||||
* is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
|
||||
*/
|
||||
public MidPointIntegrator(final int minimalIterationCount,
|
||||
final int maximalIterationCount)
|
||||
throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
|
||||
super(minimalIterationCount, maximalIterationCount);
|
||||
if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
|
||||
throw new NumberIsTooLargeException(maximalIterationCount,
|
||||
MIDPOINT_MAX_ITERATIONS_COUNT, false);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Construct a midpoint integrator with default settings.
|
||||
* (max iteration count set to {@link #MIDPOINT_MAX_ITERATIONS_COUNT})
|
||||
*/
|
||||
public MidPointIntegrator() {
|
||||
super(DEFAULT_MIN_ITERATIONS_COUNT, MIDPOINT_MAX_ITERATIONS_COUNT);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the n-th stage integral of midpoint 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
|
||||
* already computed values.</p>
|
||||
*
|
||||
* @param n the stage of 1/2 refinement, n = 0 is no refinement
|
||||
* @return the value of n-th stage integral
|
||||
* @throws TooManyEvaluationsException if the maximal number of evaluations
|
||||
* is exceeded.
|
||||
*/
|
||||
private double stage(final int n)
|
||||
throws TooManyEvaluationsException {
|
||||
|
||||
final double max = getMax();
|
||||
final double min = getMin();
|
||||
|
||||
if (n == 0) {
|
||||
final double midPoint = 0.5 * (max - min);
|
||||
s = (max - min) * computeObjectiveValue(midPoint);
|
||||
return s;
|
||||
} else {
|
||||
final long np = 1L << (n - 1); // number of new points in this stage
|
||||
double sum = 0;
|
||||
// spacing between adjacent new points
|
||||
final double spacing = (max - min) / np;
|
||||
double x = min + 0.5 * spacing; // the first new point
|
||||
for (long i = 0; i < np; i++) {
|
||||
sum += computeObjectiveValue(x);
|
||||
x += spacing;
|
||||
}
|
||||
// add the new sum to previously calculated result
|
||||
s = 0.5 * (s + sum * spacing);
|
||||
return s;
|
||||
}
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
protected double doIntegrate()
|
||||
throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {
|
||||
|
||||
double oldt = stage(0);
|
||||
iterations.incrementCount();
|
||||
while (true) {
|
||||
final int i = iterations.getCount();
|
||||
final double t = stage(i);
|
||||
if (i >= getMinimalIterationCount()) {
|
||||
final double delta = FastMath.abs(t - oldt);
|
||||
final double rLimit =
|
||||
getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5;
|
||||
if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
|
||||
return t;
|
||||
}
|
||||
}
|
||||
oldt = t;
|
||||
iterations.incrementCount();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,132 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math3.analysis.integration;
|
||||
|
||||
import org.apache.commons.math3.analysis.QuinticFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math3.analysis.function.Sin;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
||||
/**
|
||||
* Test case for midpoint 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 $Id: MidPointIntegratorTest.java 1374632 2012-08-18 18:11:11Z luc $
|
||||
*/
|
||||
public final class MidPointIntegratorTest {
|
||||
|
||||
/**
|
||||
* Test of integrator for the sine function.
|
||||
*/
|
||||
@Test
|
||||
public void testSinFunction() {
|
||||
UnivariateFunction f = new Sin();
|
||||
UnivariateIntegrator integrator = new MidPointIntegrator();
|
||||
|
||||
double min = 0;
|
||||
double max = FastMath.PI;
|
||||
double expected = 2;
|
||||
double tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
|
||||
double result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
|
||||
Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
|
||||
Assert.assertTrue(integrator.getIterations() < MidPointIntegrator.MIDPOINT_MAX_ITERATIONS_COUNT / 2);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
min = -FastMath.PI/3;
|
||||
max = 0;
|
||||
expected = -0.5;
|
||||
tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
|
||||
result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
|
||||
Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
|
||||
Assert.assertTrue(integrator.getIterations() < MidPointIntegrator.MIDPOINT_MAX_ITERATIONS_COUNT / 2);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Test of integrator for the quintic function.
|
||||
*/
|
||||
@Test
|
||||
public void testQuinticFunction() {
|
||||
UnivariateFunction f = new QuinticFunction();
|
||||
UnivariateIntegrator integrator = new MidPointIntegrator();
|
||||
|
||||
double min = 0;
|
||||
double max = 1;
|
||||
double expected = -1.0 / 48;
|
||||
double tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
|
||||
double result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
|
||||
Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
|
||||
Assert.assertTrue(integrator.getIterations() < MidPointIntegrator.MIDPOINT_MAX_ITERATIONS_COUNT / 2);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
min = 0;
|
||||
max = 0.5;
|
||||
expected = 11.0 / 768;
|
||||
tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
|
||||
result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
|
||||
Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
|
||||
Assert.assertTrue(integrator.getIterations() < MidPointIntegrator.MIDPOINT_MAX_ITERATIONS_COUNT / 2);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
min = -1;
|
||||
max = 4;
|
||||
expected = 2048 / 3.0 - 78 + 1.0 / 48;
|
||||
tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
|
||||
result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
|
||||
Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
|
||||
Assert.assertTrue(integrator.getIterations() < MidPointIntegrator.MIDPOINT_MAX_ITERATIONS_COUNT / 2);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Test of parameters for the integrator.
|
||||
*/
|
||||
@Test
|
||||
public void testParameters() {
|
||||
UnivariateFunction f = new Sin();
|
||||
|
||||
try {
|
||||
// bad interval
|
||||
new MidPointIntegrator().integrate(1000, f, 1, -1);
|
||||
Assert.fail("Expecting NumberIsTooLargeException - bad interval");
|
||||
} catch (NumberIsTooLargeException ex) {
|
||||
// expected
|
||||
}
|
||||
try {
|
||||
// bad iteration limits
|
||||
new MidPointIntegrator(5, 4);
|
||||
Assert.fail("Expecting NumberIsTooSmallException - bad iteration limits");
|
||||
} catch (NumberIsTooSmallException ex) {
|
||||
// expected
|
||||
}
|
||||
try {
|
||||
// bad iteration limits
|
||||
new MidPointIntegrator(10, 99);
|
||||
Assert.fail("Expecting NumberIsTooLargeException - bad iteration limits");
|
||||
} catch (NumberIsTooLargeException ex) {
|
||||
// expected
|
||||
}
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue