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