From 3cce9ed6c38c9a33f31c7d2d0f018bcc4221548d Mon Sep 17 00:00:00 2001 From: aherbert Date: Tue, 8 May 2018 10:46:05 +0100 Subject: [PATCH] MATH-1458: Added JUnit tests to document failure --- .../integration/SimpsonIntegratorTest.java | 266 +++++++++++++++++- 1 file changed, 265 insertions(+), 1 deletion(-) diff --git a/src/test/java/org/apache/commons/math4/analysis/integration/SimpsonIntegratorTest.java b/src/test/java/org/apache/commons/math4/analysis/integration/SimpsonIntegratorTest.java index 22947b323..53563c6e4 100644 --- a/src/test/java/org/apache/commons/math4/analysis/integration/SimpsonIntegratorTest.java +++ b/src/test/java/org/apache/commons/math4/analysis/integration/SimpsonIntegratorTest.java @@ -18,6 +18,8 @@ package org.apache.commons.math4.analysis.integration; import org.apache.commons.math4.analysis.QuinticFunction; import org.apache.commons.math4.analysis.UnivariateFunction; +import org.apache.commons.math4.analysis.function.Identity; +import org.apache.commons.math4.analysis.function.Inverse; import org.apache.commons.math4.analysis.function.Sin; import org.apache.commons.math4.analysis.integration.SimpsonIntegrator; import org.apache.commons.math4.analysis.integration.UnivariateIntegrator; @@ -114,10 +116,272 @@ public final class SimpsonIntegratorTest { } try { // bad iteration limits - new SimpsonIntegrator(10, 99); + new SimpsonIntegrator(10, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT + 1); Assert.fail("Expecting NumberIsTooLargeException - bad iteration limits"); } catch (NumberIsTooLargeException ex) { // expected } } + + // Tests for MATH-1458: + // The SimpsonIntegrator had the following bugs: + // - minimalIterationCount==1 results in no possible iteration + // - minimalIterationCount==1 computes incorrect Simpson sum (following no iteration) + // - minimalIterationCount>1 computes the first iteration sum as the Trapezoid sum + // - minimalIterationCount>1 computes the second iteration sum as the first Simpson sum + + /** + * Test iteration is possible when minimalIterationCount==1. + *
+ * MATH-1458: No iterations were performed when minimalIterationCount==1. + */ + @Test + public void testIterationIsPossibleWhenMinimalIterationCountIs1() { + UnivariateFunction f = new Sin(); + UnivariateIntegrator integrator = new SimpsonIntegrator(1, + SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT); + // The range or result is not relevant. + // This sum should not converge at 1 iteration. + // This tests iteration occurred. + integrator.integrate(1000, f, 0, 1); + // MATH-1458: No iterations were performed when minimalIterationCount==1 + Assert.assertTrue("Iteration is not above 1", + integrator.getIterations() > 1); + } + + /** + * Test convergence at iteration 1 when minimalIterationCount==1. + *
+ * MATH-1458: No iterations were performed when minimalIterationCount==1. + */ + @Test + public void testConvergenceIsPossibleAtIteration1() { + // A linear function y=x should converge immediately + UnivariateFunction f = new Identity(); + UnivariateIntegrator integrator = new SimpsonIntegrator(1, + SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT); + + double min, max, expected, result, tolerance; + + min = 0; max = 1; expected = 0.5; + tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy()); + result = integrator.integrate(1000, f, min, max); + // MATH-1458: No iterations were performed when minimalIterationCount==1 + Assert.assertTrue("Iteration is not above 0", + integrator.getIterations() > 0); + // This should converge immediately + Assert.assertEquals("Iteration", integrator.getIterations(), 1); + Assert.assertEquals("Result", expected, result, tolerance); + } + + /** + * Compute the integral using the composite Simpson's rule. + * + * @param f the function + * @param a the lower limit + * @param b the upper limit + * @param n the number of intervals (must be even) + * @return the integral between a and b + * @see + * Composite_Simpson's_rule + */ + private static double compositeSimpsonsRule(UnivariateFunction f, double a, + double b, int n) + { + // Sum interval [a,b] split into n subintervals, with n an even number: + // sum ~ h/3 * [ f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + 2f(x4) ... + 4f(xn-1) + f(xn) ] + // h = (b-a)/n + // f(xi) = f(a + i*h) + assert n > 0 && n % 2 == 0 : "n must be strictly positive and even"; + final double h = (b - a) / n; + double sum4 = 0; + double sum2 = 0; + for (int i = 1; i < n; i++) { + // Alternate sums that are multiplied by 4 and 2 + final double fxi = f.value(a + i * h); + if (i % 2 == 0) + sum2 += fxi; + else + sum4 += fxi; + } + return (h / 3) * (f.value(a) + 4 * sum4 + 2 * sum2 + f.value(b)); + } + + /** + * Compute the iteration of Simpson's rule. + * + * @param f the function + * @param a the lower limit + * @param b the upper limit + * @param iteration the refinement iteration + * @return the integral between a and b + */ + private static double computeSimpsonIteration(UnivariateFunction f, double a, + double b, int iteration) + { + // The first possible Simpson's sum uses n=2. + // The next uses n=4. This is the 1st refinement expected when the + // integrator has performed 1 iteration. + final int n = 2 << iteration; + return compositeSimpsonsRule(f, a, b, n); + } + + /** + * Test the reference Simpson integration is doing what is expected + */ + @Test + public void testReferenceSimpsonItegrationIsCorrect() { + UnivariateFunction f = new Sin(); + + double a, b, h, expected, result, tolerance; + + a = 0.5; + b = 1; + + double b_a = b - a; + + // First Simpson sum. 1 midpoint evaluation: + h = b_a / 2; + double f00 = f.value(a); + double f01 = f.value(a + 1 * h); + double f0n = f.value(b); + expected = (b_a / 6) * (f00 + 4 * f01 + f0n); + tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY); + result = computeSimpsonIteration(f, a, b, 0); + Assert.assertEquals("Result", expected, result, tolerance); + + // Second Simpson sum: 2 more evaluations: + h = b_a / 4; + double f11 = f.value(a + 1 * h); + double f13 = f.value(a + 3 * h); + expected = (h / 3) * (f00 + 4 * f11 + 2 * f01 + 4 * f13 + f0n); + tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY); + result = computeSimpsonIteration(f, a, b, 1); + Assert.assertEquals("Result", expected, result, tolerance); + + // Third Simpson sum: 4 more evaluations: + h = b_a / 8; + double f21 = f.value(a + 1 * h); + double f23 = f.value(a + 3 * h); + double f25 = f.value(a + 5 * h); + double f27 = f.value(a + 7 * h); + expected = (h / 3) * (f00 + 4 * f21 + 2 * f11 + 4 * f23 + 2 * f01 + 4 * f25 + + 2 * f13 + 4 * f27 + f0n); + tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY); + result = computeSimpsonIteration(f, a, b, 2); + Assert.assertEquals("Result", expected, result, tolerance); + } + + /** + * Test iteration 1 returns the expected sum when minimalIterationCount==1. + *
+ * MATH-1458: minimalIterationCount==1 computes incorrect Simpson sum + * (following no iteration). + */ + @Test + public void testIteration1ComputesTheExpectedSimpsonSum() { + UnivariateFunction f = new Sin(); + // Set convergence criteria to force immediate convergence + UnivariateIntegrator integrator = new SimpsonIntegrator( + 0, Double.POSITIVE_INFINITY, + 1, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT); + double min, max, expected, result, tolerance; + + // MATH-1458: minimalIterationCount==1 computes incorrect + // Simpson sum (following no iteration) + min = 0; + max = 1; + result = integrator.integrate(1000, f, min, max); + // Immediate convergence + Assert.assertEquals("Iteration", 1, integrator.getIterations()); + + // Check the sum is as expected + expected = computeSimpsonIteration(f, min, max, 1); + tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY); + Assert.assertEquals("Result", expected, result, tolerance); + } + + /** + * Test iteration N returns the expected sum when minimalIterationCount==1. + *
+ * MATH-1458: minimalIterationCount>1 computes the second iteration sum as + * the first Simpson sum. + */ + @Test + public void testIterationNComputesTheExpectedSimpsonSum() { + // Use 1/x as the function as the sum will asymptote in a monotonic + // series. The convergence can then be controlled. + UnivariateFunction f = new Inverse(); + + double min, max, expected, result, tolerance; + int minIteration, maxIteration; + + // Range for integration + min = 1; + max = 2; + + // This is the expected sum. + // Each iteration will monotonically converge to this. + expected = FastMath.log(max) - FastMath.log(min); + + // Test convergence at the given iteration + minIteration = 2; + maxIteration = 4; + + // Compute the sums expected for different iterations. + // Add an additional sum so that the test can compare to the next value. + double[] sums = new double[maxIteration + 2]; + for (int i = 0; i < sums.length; i++) { + sums[i] = computeSimpsonIteration(f, min, max, i); + // Check monotonic + if (i > 0) { + Assert.assertTrue("Expected series not monotonic descending", + sums[i] < sums[i - 1]); + // Check monotonic difference + if (i > 1) { + Assert.assertTrue("Expected convergence not monotonic descending", + sums[i - 1] - sums[i] < sums[i - 2] - sums[i - 1]); + } + } + } + + // Check the test function is correct. + tolerance = FastMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY); + Assert.assertEquals("Expected result", expected, sums[maxIteration], tolerance); + + // Set-up to test convergence at a specific iteration. + // Allow enough function evaluations. + // Iteration 0 = 3 evaluations + // Iteration 1 = 5 evaluations + // Iteration n = 2^(n+1)+1 evaluations + int evaluations = 2 << (maxIteration + 1) + 1; + + for (int i = minIteration; i <= maxIteration; i++) { + // Create convergence criteria. + // (sum - previous) is monotonic descending. + // So use a point half-way between them: + // ((sums[i-1] - sums[i]) + (sums[i-2] - sums[i-1])) / 2 + final double absoluteAccuracy = (sums[i - 2] - sums[i]) / 2; + + // Use minimalIterationCount>1 + UnivariateIntegrator integrator = new SimpsonIntegrator( + 0, absoluteAccuracy, + 2, SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT); + + result = integrator.integrate(evaluations, f, min, max); + + // Check the iteration is as expected + Assert.assertEquals("Test failed to control iteration", i, integrator.getIterations()); + + // MATH-1458: minimalIterationCount>1 computes incorrect Simpson sum + // for the iteration. Check it is the correct sum. + // It should be closer to this one than the previous or next. + final double dp = FastMath.abs(sums[i-1] - result); + final double d = FastMath.abs(sums[i] - result); + final double dn = FastMath.abs(sums[i+1] - result); + + Assert.assertTrue("Result closer to sum expected from previous iteration: " + i, d < dp); + Assert.assertTrue("Result closer to sum expected from next iteration: " + i, d < dn); + } + } }