1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.analysis.integration;
23
24 import java.util.Random;
25
26 import org.hipparchus.analysis.QuinticFunction;
27 import org.hipparchus.analysis.UnivariateFunction;
28 import org.hipparchus.analysis.function.Gaussian;
29 import org.hipparchus.analysis.function.Sin;
30 import org.hipparchus.analysis.polynomials.PolynomialFunction;
31 import org.hipparchus.exception.LocalizedCoreFormats;
32 import org.hipparchus.exception.MathIllegalStateException;
33 import org.hipparchus.util.FastMath;
34 import org.junit.Assert;
35 import org.junit.Test;
36
37
38 public class IterativeLegendreGaussIntegratorTest {
39
40 @Test
41 public void testSinFunction() {
42 UnivariateFunction f = new Sin();
43 BaseAbstractUnivariateIntegrator integrator
44 = new IterativeLegendreGaussIntegrator(5, 1.0e-14, 1.0e-10, 2, 15);
45 double min, max, expected, result, tolerance;
46
47 min = 0; max = FastMath.PI; expected = 2;
48 tolerance = FastMath.max(integrator.getAbsoluteAccuracy(),
49 FastMath.abs(expected * integrator.getRelativeAccuracy()));
50 result = integrator.integrate(10000, f, min, max);
51 Assert.assertEquals(expected, result, tolerance);
52
53 min = -FastMath.PI/3; max = 0; expected = -0.5;
54 tolerance = FastMath.max(integrator.getAbsoluteAccuracy(),
55 FastMath.abs(expected * integrator.getRelativeAccuracy()));
56 result = integrator.integrate(10000, f, min, max);
57 Assert.assertEquals(expected, result, tolerance);
58 }
59
60 @Test
61 public void testQuinticFunction() {
62 UnivariateFunction f = new QuinticFunction();
63 UnivariateIntegrator integrator =
64 new IterativeLegendreGaussIntegrator(3,
65 BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
66 BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
67 BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
68 64);
69 double min, max, expected, result;
70
71 min = 0; max = 1; expected = -1.0/48;
72 result = integrator.integrate(10000, f, min, max);
73 Assert.assertEquals(expected, result, 1.0e-16);
74
75 min = 0; max = 0.5; expected = 11.0/768;
76 result = integrator.integrate(10000, f, min, max);
77 Assert.assertEquals(expected, result, 1.0e-16);
78
79 min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
80 result = integrator.integrate(10000, f, min, max);
81 Assert.assertEquals(expected, result, 1.0e-16);
82 }
83
84 @Test
85 public void testExactIntegration() {
86 Random random = new Random(86343623467878363l);
87 for (int n = 2; n < 6; ++n) {
88 IterativeLegendreGaussIntegrator integrator =
89 new IterativeLegendreGaussIntegrator(n,
90 BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY,
91 BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY,
92 BaseAbstractUnivariateIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
93 64);
94
95
96 for (int degree = 0; degree <= 2 * n - 1; ++degree) {
97 for (int i = 0; i < 10; ++i) {
98 double[] coeff = new double[degree + 1];
99 for (int k = 0; k < coeff.length; ++k) {
100 coeff[k] = 2 * random.nextDouble() - 1;
101 }
102 PolynomialFunction p = new PolynomialFunction(coeff);
103 double result = integrator.integrate(10000, p, -5.0, 15.0);
104 double reference = exactIntegration(p, -5.0, 15.0);
105 Assert.assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + FastMath.abs(reference)));
106 }
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 * FastMath.sqrt(2 * FastMath.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(50, 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 MathIllegalStateException");
152 } catch (MathIllegalStateException tmee) {
153
154 Assert.assertEquals(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, tmee.getSpecifier());
155 Assert.assertEquals(1000, ((Integer) tmee.getParts()[0]).intValue());
156 }
157
158
159 double sum1 = gauss.integrate(1000, f, -10, 0);
160 int eval1 = gauss.getEvaluations();
161 double sum2 = gauss.integrate(1000, f, 0, maxX);
162 int eval2 = gauss.getEvaluations();
163 Assert.assertEquals(maxX * value, sum1 + sum2, 1.0e-7);
164 Assert.assertTrue(eval1 + eval2 < 200);
165
166 }
167
168 private double exactIntegration(PolynomialFunction p, double a, double b) {
169 final double[] coeffs = p.getCoefficients();
170 double yb = coeffs[coeffs.length - 1] / coeffs.length;
171 double ya = yb;
172 for (int i = coeffs.length - 2; i >= 0; --i) {
173 yb = yb * b + coeffs[i] / (i + 1);
174 ya = ya * a + coeffs[i] / (i + 1);
175 }
176 return yb * b - ya * a;
177 }
178 }