commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject 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
Date Sun, 22 Jul 2012 21:56:20 GMT
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



Mime
View raw message