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 org.apache.commons.math4.legacy.analysis.QuinticFunction;
20 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
21 import org.apache.commons.math4.legacy.analysis.function.Identity;
22 import org.apache.commons.math4.legacy.analysis.function.Inverse;
23 import org.apache.commons.math4.legacy.analysis.function.Sin;
24 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
25 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
26 import org.apache.commons.math4.core.jdkmath.JdkMath;
27 import org.junit.Assert;
28 import org.junit.Test;
29
30
31
32
33
34
35
36
37
38 public final class SimpsonIntegratorTest {
39 private static final int SIMPSON_MAX_ITERATIONS_COUNT = 30;
40
41
42
43
44 @Test
45 public void testSinFunction() {
46 UnivariateFunction f = new Sin();
47 UnivariateIntegrator integrator = new SimpsonIntegrator();
48 double min;
49 double max;
50 double expected;
51 double result;
52 double tolerance;
53
54 min = 0; max = JdkMath.PI; expected = 2;
55 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
56 result = integrator.integrate(1000, f, min, max);
57 Assert.assertTrue(integrator.getEvaluations() < 100);
58 Assert.assertTrue(integrator.getIterations() < 10);
59 Assert.assertEquals(expected, result, tolerance);
60
61 min = -JdkMath.PI/3; max = 0; expected = -0.5;
62 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
63 result = integrator.integrate(1000, f, min, max);
64 Assert.assertTrue(integrator.getEvaluations() < 50);
65 Assert.assertTrue(integrator.getIterations() < 10);
66 Assert.assertEquals(expected, result, tolerance);
67 }
68
69
70
71
72 @Test
73 public void testQuinticFunction() {
74 UnivariateFunction f = new QuinticFunction();
75 UnivariateIntegrator integrator = new SimpsonIntegrator();
76 double min;
77 double max;
78 double expected;
79 double result;
80 double tolerance;
81
82 min = 0; max = 1; expected = -1.0/48;
83 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
84 result = integrator.integrate(1000, f, min, max);
85 Assert.assertTrue(integrator.getEvaluations() < 150);
86 Assert.assertTrue(integrator.getIterations() < 10);
87 Assert.assertEquals(expected, result, tolerance);
88
89 min = 0; max = 0.5; expected = 11.0/768;
90 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
91 result = integrator.integrate(1000, f, min, max);
92 Assert.assertTrue(integrator.getEvaluations() < 100);
93 Assert.assertTrue(integrator.getIterations() < 10);
94 Assert.assertEquals(expected, result, tolerance);
95
96 min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
97 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
98 result = integrator.integrate(1000, f, min, max);
99 Assert.assertTrue(integrator.getEvaluations() < 150);
100 Assert.assertTrue(integrator.getIterations() < 10);
101 Assert.assertEquals(expected, result, tolerance);
102 }
103
104
105
106
107 @Test
108 public void testParameters() {
109 UnivariateFunction f = new Sin();
110 try {
111
112 new SimpsonIntegrator().integrate(1000, f, 1, -1);
113 Assert.fail("Expecting NumberIsTooLargeException - bad interval");
114 } catch (NumberIsTooLargeException ex) {
115
116 }
117 try {
118
119 new SimpsonIntegrator(5, 4);
120 Assert.fail("Expecting NumberIsTooSmallException - bad iteration limits");
121 } catch (NumberIsTooSmallException ex) {
122
123 }
124 try {
125
126 new SimpsonIntegrator(10, SIMPSON_MAX_ITERATIONS_COUNT + 1);
127 Assert.fail("Expecting NumberIsTooLargeException - bad iteration limits");
128 } catch (NumberIsTooLargeException ex) {
129
130 }
131 }
132
133
134
135
136
137
138
139
140
141
142
143
144
145 @Test
146 public void testIterationIsPossibleWhenMinimalIterationCountIs1() {
147 UnivariateFunction f = new Sin();
148 UnivariateIntegrator integrator = new SimpsonIntegrator(1, SIMPSON_MAX_ITERATIONS_COUNT);
149
150
151
152 integrator.integrate(1000, f, 0, 1);
153
154 Assert.assertTrue("Iteration is not above 1",
155 integrator.getIterations() > 1);
156 }
157
158
159
160
161
162
163 @Test
164 public void testConvergenceIsPossibleAtIteration1() {
165
166 UnivariateFunction f = new Identity();
167 UnivariateIntegrator integrator = new SimpsonIntegrator(1, SIMPSON_MAX_ITERATIONS_COUNT);
168
169 double min;
170 double max;
171 double expected;
172 double result;
173 double tolerance;
174
175 min = 0; max = 1; expected = 0.5;
176 tolerance = JdkMath.abs(expected * integrator.getRelativeAccuracy());
177 result = integrator.integrate(1000, f, min, max);
178
179 Assert.assertTrue("Iteration is not above 0",
180 integrator.getIterations() > 0);
181
182 Assert.assertEquals("Iteration", integrator.getIterations(), 1);
183 Assert.assertEquals("Result", expected, result, tolerance);
184 }
185
186
187
188
189
190
191
192
193
194
195
196
197 private static double compositeSimpsonsRule(UnivariateFunction f, double a,
198 double b, int n) {
199
200
201
202
203 assert n > 0 && (n & 1) == 0 : "n must be strictly positive and even";
204 final double h = (b - a) / n;
205 double sum4 = 0;
206 double sum2 = 0;
207 for (int i = 1; i < n; i++) {
208
209 final double fxi = f.value(a + i * h);
210 if ((i & 1) == 0) {
211 sum2 += fxi;
212 } else {
213 sum4 += fxi;
214 }
215 }
216 return (h / 3) * (f.value(a) + 4 * sum4 + 2 * sum2 + f.value(b));
217 }
218
219
220
221
222
223
224
225
226
227
228 private static double computeSimpsonIteration(UnivariateFunction f, double a,
229 double b, int iteration) {
230
231
232
233 final int n = 2 << iteration;
234 return compositeSimpsonsRule(f, a, b, n);
235 }
236
237
238
239
240 @Test
241 public void testReferenceSimpsonItegrationIsCorrect() {
242 UnivariateFunction f = new Sin();
243
244 double a;
245 double b;
246 double h;
247 double expected;
248 double result;
249 double tolerance;
250
251 a = 0.5;
252 b = 1;
253
254 double b_a = b - a;
255
256
257 h = b_a / 2;
258 double f00 = f.value(a);
259 double f01 = f.value(a + 1 * h);
260 double f0n = f.value(b);
261 expected = (b_a / 6) * (f00 + 4 * f01 + f0n);
262 tolerance = JdkMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
263 result = computeSimpsonIteration(f, a, b, 0);
264 Assert.assertEquals("Result", expected, result, tolerance);
265
266
267 h = b_a / 4;
268 double f11 = f.value(a + 1 * h);
269 double f13 = f.value(a + 3 * h);
270 expected = (h / 3) * (f00 + 4 * f11 + 2 * f01 + 4 * f13 + f0n);
271 tolerance = JdkMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
272 result = computeSimpsonIteration(f, a, b, 1);
273 Assert.assertEquals("Result", expected, result, tolerance);
274
275
276 h = b_a / 8;
277 double f21 = f.value(a + 1 * h);
278 double f23 = f.value(a + 3 * h);
279 double f25 = f.value(a + 5 * h);
280 double f27 = f.value(a + 7 * h);
281 expected = (h / 3) * (f00 + 4 * f21 + 2 * f11 + 4 * f23 + 2 * f01 + 4 * f25 +
282 2 * f13 + 4 * f27 + f0n);
283 tolerance = JdkMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
284 result = computeSimpsonIteration(f, a, b, 2);
285 Assert.assertEquals("Result", expected, result, tolerance);
286 }
287
288
289
290
291
292
293
294 @Test
295 public void testIteration1ComputesTheExpectedSimpsonSum() {
296 UnivariateFunction f = new Sin();
297
298 UnivariateIntegrator integrator = new SimpsonIntegrator(
299 0, Double.POSITIVE_INFINITY,
300 1, SIMPSON_MAX_ITERATIONS_COUNT);
301 double min;
302 double max;
303 double expected;
304 double result;
305 double tolerance;
306
307
308
309 min = 0;
310 max = 1;
311 result = integrator.integrate(1000, f, min, max);
312
313 Assert.assertEquals("Iteration", 1, integrator.getIterations());
314
315
316 expected = computeSimpsonIteration(f, min, max, 1);
317 tolerance = JdkMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
318 Assert.assertEquals("Result", expected, result, tolerance);
319 }
320
321
322
323
324
325
326
327 @Test
328 public void testIterationNComputesTheExpectedSimpsonSum() {
329
330
331 UnivariateFunction f = new Inverse();
332
333 double min;
334 double max;
335 double expected;
336 double result;
337 double tolerance;
338 int minIteration;
339 int maxIteration;
340
341
342 min = 1;
343 max = 2;
344
345
346
347 expected = JdkMath.log(max) - JdkMath.log(min);
348
349
350 minIteration = 2;
351 maxIteration = 4;
352
353
354
355 double[] sums = new double[maxIteration + 2];
356 for (int i = 0; i < sums.length; i++) {
357 sums[i] = computeSimpsonIteration(f, min, max, i);
358
359 if (i > 0) {
360 Assert.assertTrue("Expected series not monotonic descending",
361 sums[i] < sums[i - 1]);
362
363 if (i > 1) {
364 Assert.assertTrue("Expected convergence not monotonic descending",
365 sums[i - 1] - sums[i] < sums[i - 2] - sums[i - 1]);
366 }
367 }
368 }
369
370
371 tolerance = JdkMath.abs(expected * SimpsonIntegrator.DEFAULT_RELATIVE_ACCURACY);
372 Assert.assertEquals("Expected result", expected, sums[maxIteration], tolerance);
373
374
375
376
377
378
379 int evaluations = 2 << (maxIteration + 1) + 1;
380
381 for (int i = minIteration; i <= maxIteration; i++) {
382
383
384
385
386 final double absoluteAccuracy = (sums[i - 2] - sums[i]) / 2;
387
388
389 UnivariateIntegrator integrator = new SimpsonIntegrator(
390 0, absoluteAccuracy,
391 2, SIMPSON_MAX_ITERATIONS_COUNT);
392
393 result = integrator.integrate(evaluations, f, min, max);
394
395
396 Assert.assertEquals("Test failed to control iteration", i, integrator.getIterations());
397
398
399
400
401 final double dp = JdkMath.abs(sums[i-1] - result);
402 final double d = JdkMath.abs(sums[i] - result);
403 final double dn = JdkMath.abs(sums[i+1] - result);
404
405 Assert.assertTrue("Result closer to sum expected from previous iteration: " + i, d < dp);
406 Assert.assertTrue("Result closer to sum expected from next iteration: " + i, d < dn);
407 }
408 }
409 }