1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.analysis.integration;
18
19 import java.util.Random;
20
21 import org.apache.commons.math4.legacy.analysis.QuinticFunction;
22 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
23 import org.apache.commons.math4.legacy.analysis.function.Gaussian;
24 import org.apache.commons.math4.legacy.analysis.function.Sin;
25 import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
26 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
27 import org.apache.commons.math4.core.jdkmath.JdkMath;
28 import org.junit.Assert;
29 import org.junit.Test;
30
31
32 public class IterativeLegendreGaussIntegratorTest {
33
34 @Test
35 public void testSinFunction() {
36 UnivariateFunction f = new Sin();
37 BaseAbstractUnivariateIntegrator integrator
38 = new IterativeLegendreGaussIntegrator(5, 1.0e-14, 1.0e-10, 2, 15);
39 double min;
40 double max;
41 double expected;
42 double result;
43 double tolerance;
44
45 min = 0; max = JdkMath.PI; expected = 2;
46 tolerance = JdkMath.max(integrator.getAbsoluteAccuracy(),
47 JdkMath.abs(expected * integrator.getRelativeAccuracy()));
48 result = integrator.integrate(10000, f, min, max);
49 Assert.assertEquals(expected, result, tolerance);
50
51 min = -JdkMath.PI/3; max = 0; expected = -0.5;
52 tolerance = JdkMath.max(integrator.getAbsoluteAccuracy(),
53 JdkMath.abs(expected * integrator.getRelativeAccuracy()));
54 result = integrator.integrate(10000, f, min, max);
55 Assert.assertEquals(expected, result, tolerance);
56 }
57
58 @Test
59 public void testQuinticFunction() {
60 UnivariateFunction f = new QuinticFunction();
61 UnivariateIntegrator integrator =
62 new IterativeLegendreGaussIntegrator(3,
63 BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
64 BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
65 BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
66 64);
67 double min;
68 double max;
69 double expected;
70 double result;
71
72 min = 0; max = 1; expected = -1.0/48;
73 result = integrator.integrate(10000, f, min, max);
74 Assert.assertEquals(expected, result, 1.0e-16);
75
76 min = 0; max = 0.5; expected = 11.0/768;
77 result = integrator.integrate(10000, f, min, max);
78 Assert.assertEquals(expected, result, 1.0e-16);
79
80 min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
81 result = integrator.integrate(10000, f, min, max);
82 Assert.assertEquals(expected, result, 1.0e-16);
83 }
84
85 @Test
86 public void testExactIntegration() {
87 Random random = new Random(86343623467878363L);
88 for (int n = 2; n < 6; ++n) {
89 IterativeLegendreGaussIntegrator integrator =
90 new IterativeLegendreGaussIntegrator(n,
91 BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
92 BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
93 BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
94 64);
95
96
97 for (int degree = 0; degree <= 2 * n - 1; ++degree) {
98 for (int i = 0; i < 10; ++i) {
99 double[] coeff = new double[degree + 1];
100 for (int k = 0; k < coeff.length; ++k) {
101 coeff[k] = 2 * random.nextDouble() - 1;
102 }
103 PolynomialFunction p = new PolynomialFunction(coeff);
104 double result = integrator.integrate(10000, p, -5.0, 15.0);
105 double reference = exactIntegration(p, -5.0, 15.0);
106 Assert.assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + JdkMath.abs(reference)));
107 }
108 }
109 }
110 }
111
112
113 @Test
114 public void testNormalDistributionWithLargeSigma() {
115 final double sigma = 1000;
116 final double mean = 0;
117 final double factor = 1 / (sigma * JdkMath.sqrt(2 * JdkMath.PI));
118 final UnivariateFunction normal = new Gaussian(factor, mean, sigma);
119
120 final double tol = 1e-2;
121 final IterativeLegendreGaussIntegrator integrator =
122 new IterativeLegendreGaussIntegrator(5, tol, tol);
123
124 final double a = -5000;
125 final double b = 5000;
126 final double s = integrator.integrate(60, normal, a, b);
127 Assert.assertEquals(1, s, 1e-5);
128 }
129
130 @Test
131 public void testIssue464() {
132 final double value = 0.2;
133 UnivariateFunction f = new UnivariateFunction() {
134 @Override
135 public double value(double x) {
136 return (x >= 0 && x <= 5) ? value : 0.0;
137 }
138 };
139 IterativeLegendreGaussIntegrator gauss
140 = new IterativeLegendreGaussIntegrator(5, 3, 100);
141
142
143 double maxX = 0.32462367623786328;
144 Assert.assertEquals(maxX * value, gauss.integrate(Integer.MAX_VALUE, f, -10, maxX), 1.0e-7);
145 Assert.assertTrue(gauss.getEvaluations() > 37000000);
146 Assert.assertTrue(gauss.getIterations() < 30);
147
148
149 try {
150 gauss.integrate(1000, f, -10, maxX);
151 Assert.fail("expected TooManyEvaluationsException");
152 } catch (TooManyEvaluationsException tmee) {
153
154 Assert.assertEquals(1000, tmee.getMax());
155 }
156
157
158 double sum1 = gauss.integrate(1000, f, -10, 0);
159 int eval1 = gauss.getEvaluations();
160 double sum2 = gauss.integrate(1000, f, 0, maxX);
161 int eval2 = gauss.getEvaluations();
162 Assert.assertEquals(maxX * value, sum1 + sum2, 1.0e-7);
163 Assert.assertTrue(eval1 + eval2 < 200);
164 }
165
166 private double exactIntegration(PolynomialFunction p, double a, double b) {
167 final double[] coeffs = p.getCoefficients();
168 double yb = coeffs[coeffs.length - 1] / coeffs.length;
169 double ya = yb;
170 for (int i = coeffs.length - 2; i >= 0; --i) {
171 yb = yb * b + coeffs[i] / (i + 1);
172 ya = ya * a + coeffs[i] / (i + 1);
173 }
174 return yb * b - ya * a;
175 }
176 }