MATH-827
New "IterativeLegendreGaussIntegrator" class that performs the same algorithm as the current "LegendreGaussIntegrator" but uses the recently added Gauss integration framework (in package "o.a.c.m.analysis.integration.gauss") for the underlying integration computations. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1364444 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski
parent
463896f1cc
commit
28a66efb76
|
|
@ -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.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory;
|
||||
import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator;
|
||||
import org.apache.commons.math3.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math3.exception.TooManyEvaluationsException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
|
||||
/**
|
||||
* This algorithm divides the integration interval into equally-sized
|
||||
* sub-interval and on each of them performs a
|
||||
* <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
|
||||
* Legendre-Gauss</a> quadrature.
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 3.1
|
||||
*/
|
||||
|
||||
public class IterativeLegendreGaussIntegrator
|
||||
extends BaseAbstractUnivariateIntegrator {
|
||||
/** Factory that computes the points and weights. */
|
||||
private static final GaussIntegratorFactory factory
|
||||
= new GaussIntegratorFactory();
|
||||
/** Number of integration points (per interval). */
|
||||
private final int numberOfPoints;
|
||||
|
||||
/**
|
||||
* Builds an integrator with given accuracies and iterations counts.
|
||||
*
|
||||
* @param n Number of integration points.
|
||||
* @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.
|
||||
* @throws NotStrictlyPositiveException if minimal number of iterations
|
||||
* is not strictly positive.
|
||||
* @throws NumberIsTooSmallException if maximal number of iterations
|
||||
* is smaller than or equal to the minimal number of iterations.
|
||||
*/
|
||||
public IterativeLegendreGaussIntegrator(final int n,
|
||||
final double relativeAccuracy,
|
||||
final double absoluteAccuracy,
|
||||
final int minimalIterationCount,
|
||||
final int maximalIterationCount)
|
||||
throws NotStrictlyPositiveException, NumberIsTooSmallException {
|
||||
super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
|
||||
numberOfPoints = n;
|
||||
}
|
||||
|
||||
/**
|
||||
* Builds an integrator with given accuracies.
|
||||
*
|
||||
* @param n Number of integration points.
|
||||
* @param relativeAccuracy Relative accuracy of the result.
|
||||
* @param absoluteAccuracy Absolute accuracy of the result.
|
||||
*/
|
||||
public IterativeLegendreGaussIntegrator(final int n,
|
||||
final double relativeAccuracy,
|
||||
final double absoluteAccuracy) {
|
||||
this(n, relativeAccuracy, absoluteAccuracy,
|
||||
DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
|
||||
}
|
||||
|
||||
/**
|
||||
* Builds an integrator with given iteration counts.
|
||||
*
|
||||
* @param n Number of integration points.
|
||||
* @param minimalIterationCount Minimum number of iterations.
|
||||
* @param maximalIterationCount Maximum number of iterations.
|
||||
* @throws NotStrictlyPositiveException if minimal number of iterations
|
||||
* is not strictly positive.
|
||||
* @throws NumberIsTooSmallException if maximal number of iterations
|
||||
* is smaller than or equal to the minimal number of iterations.
|
||||
*/
|
||||
public IterativeLegendreGaussIntegrator(final int n,
|
||||
final int minimalIterationCount,
|
||||
final int maximalIterationCount) {
|
||||
this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
|
||||
minimalIterationCount, maximalIterationCount);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected double doIntegrate()
|
||||
throws TooManyEvaluationsException, MaxCountExceededException {
|
||||
// Compute first estimate with a single step.
|
||||
double oldt = stage(1);
|
||||
|
||||
int n = 2;
|
||||
while (true) {
|
||||
// Improve integral with a larger number of steps.
|
||||
final double t = stage(n);
|
||||
|
||||
// Estimate the error.
|
||||
final double delta = FastMath.abs(t - oldt);
|
||||
final double limit =
|
||||
FastMath.max(getAbsoluteAccuracy(),
|
||||
getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
|
||||
|
||||
// check convergence
|
||||
if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
|
||||
delta <= limit) {
|
||||
return t;
|
||||
}
|
||||
|
||||
// Prepare next iteration.
|
||||
final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
|
||||
n = FastMath.max((int) (ratio * n), n + 1);
|
||||
oldt = t;
|
||||
iterations.incrementCount();
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the n-th stage integral.
|
||||
*
|
||||
* @param n Number of steps.
|
||||
* @return the value of n-th stage integral.
|
||||
* @throws TooManyEvaluationsException if the maximum number of evaluations
|
||||
* is exceeded.
|
||||
*/
|
||||
private double stage(final int n)
|
||||
throws TooManyEvaluationsException {
|
||||
// Function to be integrated is stored in the base class.
|
||||
final UnivariateFunction f = new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
return computeObjectiveValue(x);
|
||||
}
|
||||
};
|
||||
|
||||
final double min = getMin();
|
||||
final double max = getMax();
|
||||
final double step = (max - min) / n;
|
||||
|
||||
double sum = 0;
|
||||
for (int i = 0; i < n; i++) {
|
||||
// Integrate over each sub-interval [a, b].
|
||||
final double a = min + i * step;
|
||||
final double b = a + step;
|
||||
final GaussIntegrator g = factory.legendreHighPrecision(numberOfPoints, a, b);
|
||||
sum += g.integrate(f);
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,151 @@
|
|||
/*
|
||||
* 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 java.util.Random;
|
||||
|
||||
import org.apache.commons.math3.analysis.QuinticFunction;
|
||||
import org.apache.commons.math3.analysis.SinFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
|
||||
import org.apache.commons.math3.exception.TooManyEvaluationsException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
||||
|
||||
public class IterativeLegendreGaussIntegratorTest {
|
||||
|
||||
@Test
|
||||
public void testSinFunction() {
|
||||
UnivariateFunction f = new SinFunction();
|
||||
BaseAbstractUnivariateIntegrator integrator
|
||||
= new IterativeLegendreGaussIntegrator(5, 1.0e-14, 1.0e-10, 2, 15);
|
||||
double min, max, expected, result, tolerance;
|
||||
|
||||
min = 0; max = FastMath.PI; expected = 2;
|
||||
tolerance = FastMath.max(integrator.getAbsoluteAccuracy(),
|
||||
FastMath.abs(expected * integrator.getRelativeAccuracy()));
|
||||
result = integrator.integrate(10000, f, min, max);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
|
||||
min = -FastMath.PI/3; max = 0; expected = -0.5;
|
||||
tolerance = FastMath.max(integrator.getAbsoluteAccuracy(),
|
||||
FastMath.abs(expected * integrator.getRelativeAccuracy()));
|
||||
result = integrator.integrate(10000, f, min, max);
|
||||
Assert.assertEquals(expected, result, tolerance);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testQuinticFunction() {
|
||||
UnivariateFunction f = new QuinticFunction();
|
||||
UnivariateIntegrator integrator =
|
||||
new IterativeLegendreGaussIntegrator(3,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
|
||||
64);
|
||||
double min, max, expected, result;
|
||||
|
||||
min = 0; max = 1; expected = -1.0/48;
|
||||
result = integrator.integrate(10000, f, min, max);
|
||||
Assert.assertEquals(expected, result, 1.0e-16);
|
||||
|
||||
min = 0; max = 0.5; expected = 11.0/768;
|
||||
result = integrator.integrate(10000, f, min, max);
|
||||
Assert.assertEquals(expected, result, 1.0e-16);
|
||||
|
||||
min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
|
||||
result = integrator.integrate(10000, f, min, max);
|
||||
Assert.assertEquals(expected, result, 1.0e-16);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testExactIntegration() {
|
||||
Random random = new Random(86343623467878363l);
|
||||
for (int n = 2; n < 6; ++n) {
|
||||
IterativeLegendreGaussIntegrator integrator =
|
||||
new IterativeLegendreGaussIntegrator(n,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
|
||||
BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
|
||||
64);
|
||||
|
||||
// an n points Gauss-Legendre integrator integrates 2n-1 degree polynoms exactly
|
||||
for (int degree = 0; degree <= 2 * n - 1; ++degree) {
|
||||
for (int i = 0; i < 10; ++i) {
|
||||
double[] coeff = new double[degree + 1];
|
||||
for (int k = 0; k < coeff.length; ++k) {
|
||||
coeff[k] = 2 * random.nextDouble() - 1;
|
||||
}
|
||||
PolynomialFunction p = new PolynomialFunction(coeff);
|
||||
double result = integrator.integrate(10000, p, -5.0, 15.0);
|
||||
double reference = exactIntegration(p, -5.0, 15.0);
|
||||
Assert.assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + FastMath.abs(reference)));
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testIssue464() {
|
||||
final double value = 0.2;
|
||||
UnivariateFunction f = new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
return (x >= 0 && x <= 5) ? value : 0.0;
|
||||
}
|
||||
};
|
||||
IterativeLegendreGaussIntegrator gauss
|
||||
= new IterativeLegendreGaussIntegrator(5, 3, 100);
|
||||
|
||||
// due to the discontinuity, integration implies *many* calls
|
||||
double maxX = 0.32462367623786328;
|
||||
Assert.assertEquals(maxX * value, gauss.integrate(Integer.MAX_VALUE, f, -10, maxX), 1.0e-7);
|
||||
Assert.assertTrue(gauss.getEvaluations() > 37000000);
|
||||
Assert.assertTrue(gauss.getIterations() < 30);
|
||||
|
||||
// setting up limits prevents such large number of calls
|
||||
try {
|
||||
gauss.integrate(1000, f, -10, maxX);
|
||||
Assert.fail("expected TooManyEvaluationsException");
|
||||
} catch (TooManyEvaluationsException tmee) {
|
||||
// expected
|
||||
Assert.assertEquals(1000, tmee.getMax());
|
||||
}
|
||||
|
||||
// integrating on the two sides should be simpler
|
||||
double sum1 = gauss.integrate(1000, f, -10, 0);
|
||||
int eval1 = gauss.getEvaluations();
|
||||
double sum2 = gauss.integrate(1000, f, 0, maxX);
|
||||
int eval2 = gauss.getEvaluations();
|
||||
Assert.assertEquals(maxX * value, sum1 + sum2, 1.0e-7);
|
||||
Assert.assertTrue(eval1 + eval2 < 200);
|
||||
|
||||
}
|
||||
|
||||
private double exactIntegration(PolynomialFunction p, double a, double b) {
|
||||
final double[] coeffs = p.getCoefficients();
|
||||
double yb = coeffs[coeffs.length - 1] / coeffs.length;
|
||||
double ya = yb;
|
||||
for (int i = coeffs.length - 2; i >= 0; --i) {
|
||||
yb = yb * b + coeffs[i] / (i + 1);
|
||||
ya = ya * a + coeffs[i] / (i + 1);
|
||||
}
|
||||
return yb * b - ya * a;
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue