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