| /* |
| * 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.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.exception.NumberIsTooLargeException; |
| import org.apache.commons.math4.exception.NumberIsTooSmallException; |
| import org.apache.commons.math4.util.FastMath; |
| import org.junit.Assert; |
| import org.junit.Test; |
| |
| |
| /** |
| * Test case for Simpson integrator. |
| * <p> |
| * Test runs show that for a default relative accuracy of 1E-6, it |
| * generally takes 5 to 10 iterations for the integral to converge. |
| * |
| */ |
| public final class SimpsonIntegratorTest { |
| |
| /** |
| * Test of integrator for the sine function. |
| */ |
| @Test |
| public void testSinFunction() { |
| UnivariateFunction f = new Sin(); |
| UnivariateIntegrator integrator = new SimpsonIntegrator(); |
| double min, max, expected, result, tolerance; |
| |
| min = 0; max = FastMath.PI; expected = 2; |
| tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy()); |
| result = integrator.integrate(1000, f, min, max); |
| Assert.assertTrue(integrator.getEvaluations() < 100); |
| Assert.assertTrue(integrator.getIterations() < 10); |
| Assert.assertEquals(expected, result, tolerance); |
| |
| min = -FastMath.PI/3; max = 0; expected = -0.5; |
| tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy()); |
| result = integrator.integrate(1000, f, min, max); |
| Assert.assertTrue(integrator.getEvaluations() < 50); |
| Assert.assertTrue(integrator.getIterations() < 10); |
| Assert.assertEquals(expected, result, tolerance); |
| } |
| |
| /** |
| * Test of integrator for the quintic function. |
| */ |
| @Test |
| public void testQuinticFunction() { |
| UnivariateFunction f = new QuinticFunction(); |
| UnivariateIntegrator integrator = new SimpsonIntegrator(); |
| double min, max, expected, result, tolerance; |
| |
| min = 0; max = 1; expected = -1.0/48; |
| tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy()); |
| result = integrator.integrate(1000, f, min, max); |
| Assert.assertTrue(integrator.getEvaluations() < 150); |
| Assert.assertTrue(integrator.getIterations() < 10); |
| Assert.assertEquals(expected, result, tolerance); |
| |
| min = 0; max = 0.5; expected = 11.0/768; |
| tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy()); |
| result = integrator.integrate(1000, f, min, max); |
| Assert.assertTrue(integrator.getEvaluations() < 100); |
| Assert.assertTrue(integrator.getIterations() < 10); |
| 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(1000, f, min, max); |
| Assert.assertTrue(integrator.getEvaluations() < 150); |
| Assert.assertTrue(integrator.getIterations() < 10); |
| Assert.assertEquals(expected, result, tolerance); |
| } |
| |
| /** |
| * Test of parameters for the integrator. |
| */ |
| @Test |
| public void testParameters() { |
| UnivariateFunction f = new Sin(); |
| try { |
| // bad interval |
| new SimpsonIntegrator().integrate(1000, f, 1, -1); |
| Assert.fail("Expecting NumberIsTooLargeException - bad interval"); |
| } catch (NumberIsTooLargeException ex) { |
| // expected |
| } |
| try { |
| // bad iteration limits |
| new SimpsonIntegrator(5, 4); |
| Assert.fail("Expecting NumberIsTooSmallException - bad iteration limits"); |
| } catch (NumberIsTooSmallException ex) { |
| // expected |
| } |
| try { |
| // bad iteration limits |
| 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. |
| * <br/> |
| * 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. |
| * <br/> |
| * 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 <a href="https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule"> |
| * Composite_Simpson's_rule</a> |
| */ |
| 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. |
| * <br/> |
| * 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. |
| * <br/> |
| * 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); |
| } |
| } |
| } |