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