You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2012/07/22 23:56:20 UTC
svn commit: r1364444 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java
test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java
Author: erans
Date: Sun Jul 22 21:56:20 2012
New Revision: 1364444
URL: http://svn.apache.org/viewvc?rev=1364444&view=rev
Log:
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.
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java (with props)
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java?rev=1364444&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java Sun Jul 22 21:56:20 2012
@@ -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;
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java?rev=1364444&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java Sun Jul 22 21:56:20 2012
@@ -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;
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegratorTest.java
------------------------------------------------------------------------------
svn:eol-style = native