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.polynomials;
23
24 import org.hipparchus.analysis.UnivariateFunction;
25 import org.hipparchus.analysis.integration.IterativeLegendreGaussIntegrator;
26 import org.hipparchus.util.CombinatoricsUtils;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.Precision;
29 import org.junit.jupiter.api.Test;
30
31 import static org.junit.jupiter.api.Assertions.assertEquals;
32 import static org.junit.jupiter.api.Assertions.assertNotEquals;
33 import static org.junit.jupiter.api.Assertions.assertTrue;
34
35
36
37
38
39 class PolynomialsUtilsTest {
40
41 @Test
42 void testFirstChebyshevPolynomials() {
43 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(3), "-3 x + 4 x^3");
44 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(2), "-1 + 2 x^2");
45 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(1), "x");
46 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(0), "1");
47
48 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(7), "-7 x + 56 x^3 - 112 x^5 + 64 x^7");
49 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(6), "-1 + 18 x^2 - 48 x^4 + 32 x^6");
50 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(5), "5 x - 20 x^3 + 16 x^5");
51 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(4), "1 - 8 x^2 + 8 x^4");
52
53 }
54
55 @Test
56 void testChebyshevBounds() {
57 for (int k = 0; k < 12; ++k) {
58 PolynomialFunction Tk = PolynomialsUtils.createChebyshevPolynomial(k);
59 for (double x = -1; x <= 1; x += 0.02) {
60 assertTrue(FastMath.abs(Tk.value(x)) < (1 + 1e-12), k + " " + Tk.value(x));
61 }
62 }
63 }
64
65 @Test
66 void testChebyshevDifferentials() {
67 for (int k = 0; k < 12; ++k) {
68
69 PolynomialFunction Tk0 = PolynomialsUtils.createChebyshevPolynomial(k);
70 PolynomialFunction Tk1 = Tk0.polynomialDerivative();
71 PolynomialFunction Tk2 = Tk1.polynomialDerivative();
72
73 PolynomialFunction g0 = new PolynomialFunction(new double[] { k * k });
74 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -1});
75 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
76
77 PolynomialFunction Tk0g0 = Tk0.multiply(g0);
78 PolynomialFunction Tk1g1 = Tk1.multiply(g1);
79 PolynomialFunction Tk2g2 = Tk2.multiply(g2);
80
81 checkNullPolynomial(Tk0g0.add(Tk1g1.add(Tk2g2)));
82
83 }
84 }
85
86 @Test
87 void testChebyshevOrthogonality() {
88 UnivariateFunction weight = new UnivariateFunction() {
89 @Override
90 public double value(double x) {
91 return 1 / FastMath.sqrt(1 - x * x);
92 }
93 };
94 for (int i = 0; i < 10; ++i) {
95 PolynomialFunction pi = PolynomialsUtils.createChebyshevPolynomial(i);
96 for (int j = 0; j <= i; ++j) {
97 PolynomialFunction pj = PolynomialsUtils.createChebyshevPolynomial(j);
98 checkOrthogonality(pi, pj, weight, -0.9999, 0.9999, 1.5, 0.03);
99 }
100 }
101 }
102
103 @Test
104 void testFirstHermitePolynomials() {
105 checkPolynomial(PolynomialsUtils.createHermitePolynomial(3), "-12 x + 8 x^3");
106 checkPolynomial(PolynomialsUtils.createHermitePolynomial(2), "-2 + 4 x^2");
107 checkPolynomial(PolynomialsUtils.createHermitePolynomial(1), "2 x");
108 checkPolynomial(PolynomialsUtils.createHermitePolynomial(0), "1");
109
110 checkPolynomial(PolynomialsUtils.createHermitePolynomial(7), "-1680 x + 3360 x^3 - 1344 x^5 + 128 x^7");
111 checkPolynomial(PolynomialsUtils.createHermitePolynomial(6), "-120 + 720 x^2 - 480 x^4 + 64 x^6");
112 checkPolynomial(PolynomialsUtils.createHermitePolynomial(5), "120 x - 160 x^3 + 32 x^5");
113 checkPolynomial(PolynomialsUtils.createHermitePolynomial(4), "12 - 48 x^2 + 16 x^4");
114
115 }
116
117 @Test
118 void testHermiteDifferentials() {
119 for (int k = 0; k < 12; ++k) {
120
121 PolynomialFunction Hk0 = PolynomialsUtils.createHermitePolynomial(k);
122 PolynomialFunction Hk1 = Hk0.polynomialDerivative();
123 PolynomialFunction Hk2 = Hk1.polynomialDerivative();
124
125 PolynomialFunction g0 = new PolynomialFunction(new double[] { 2 * k });
126 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
127 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1 });
128
129 PolynomialFunction Hk0g0 = Hk0.multiply(g0);
130 PolynomialFunction Hk1g1 = Hk1.multiply(g1);
131 PolynomialFunction Hk2g2 = Hk2.multiply(g2);
132
133 checkNullPolynomial(Hk0g0.add(Hk1g1.add(Hk2g2)));
134
135 }
136 }
137
138 @Test
139 void testHermiteOrthogonality() {
140 UnivariateFunction weight = new UnivariateFunction() {
141 @Override
142 public double value(double x) {
143 return FastMath.exp(-x * x);
144 }
145 };
146 for (int i = 0; i < 10; ++i) {
147 PolynomialFunction pi = PolynomialsUtils.createHermitePolynomial(i);
148 for (int j = 0; j <= i; ++j) {
149 PolynomialFunction pj = PolynomialsUtils.createHermitePolynomial(j);
150 checkOrthogonality(pi, pj, weight, -50, 50, 1.5, 1.0e-8);
151 }
152 }
153 }
154
155 @Test
156 void testFirstLaguerrePolynomials() {
157 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(3), 6l, "6 - 18 x + 9 x^2 - x^3");
158 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(2), 2l, "2 - 4 x + x^2");
159 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(1), 1l, "1 - x");
160 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(0), 1l, "1");
161
162 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(7), 5040l,
163 "5040 - 35280 x + 52920 x^2 - 29400 x^3"
164 + " + 7350 x^4 - 882 x^5 + 49 x^6 - x^7");
165 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(6), 720l,
166 "720 - 4320 x + 5400 x^2 - 2400 x^3 + 450 x^4"
167 + " - 36 x^5 + x^6");
168 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(5), 120l,
169 "120 - 600 x + 600 x^2 - 200 x^3 + 25 x^4 - x^5");
170 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(4), 24l,
171 "24 - 96 x + 72 x^2 - 16 x^3 + x^4");
172
173 }
174
175 @Test
176 void testLaguerreDifferentials() {
177 for (int k = 0; k < 12; ++k) {
178
179 PolynomialFunction Lk0 = PolynomialsUtils.createLaguerrePolynomial(k);
180 PolynomialFunction Lk1 = Lk0.polynomialDerivative();
181 PolynomialFunction Lk2 = Lk1.polynomialDerivative();
182
183 PolynomialFunction g0 = new PolynomialFunction(new double[] { k });
184 PolynomialFunction g1 = new PolynomialFunction(new double[] { 1, -1 });
185 PolynomialFunction g2 = new PolynomialFunction(new double[] { 0, 1 });
186
187 PolynomialFunction Lk0g0 = Lk0.multiply(g0);
188 PolynomialFunction Lk1g1 = Lk1.multiply(g1);
189 PolynomialFunction Lk2g2 = Lk2.multiply(g2);
190
191 checkNullPolynomial(Lk0g0.add(Lk1g1.add(Lk2g2)));
192
193 }
194 }
195
196 @Test
197 void testLaguerreOrthogonality() {
198 UnivariateFunction weight = new UnivariateFunction() {
199 @Override
200 public double value(double x) {
201 return FastMath.exp(-x);
202 }
203 };
204 for (int i = 0; i < 10; ++i) {
205 PolynomialFunction pi = PolynomialsUtils.createLaguerrePolynomial(i);
206 for (int j = 0; j <= i; ++j) {
207 PolynomialFunction pj = PolynomialsUtils.createLaguerrePolynomial(j);
208 checkOrthogonality(pi, pj, weight, 0.0, 100.0, 0.99999, 1.0e-13);
209 }
210 }
211 }
212
213 @Test
214 void testFirstLegendrePolynomials() {
215 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(3), 2l, "-3 x + 5 x^3");
216 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(2), 2l, "-1 + 3 x^2");
217 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(1), 1l, "x");
218 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(0), 1l, "1");
219
220 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(7), 16l, "-35 x + 315 x^3 - 693 x^5 + 429 x^7");
221 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(6), 16l, "-5 + 105 x^2 - 315 x^4 + 231 x^6");
222 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(5), 8l, "15 x - 70 x^3 + 63 x^5");
223 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(4), 8l, "3 - 30 x^2 + 35 x^4");
224
225 }
226
227 @Test
228 void testLegendreDifferentials() {
229 for (int k = 0; k < 12; ++k) {
230
231 PolynomialFunction Pk0 = PolynomialsUtils.createLegendrePolynomial(k);
232 PolynomialFunction Pk1 = Pk0.polynomialDerivative();
233 PolynomialFunction Pk2 = Pk1.polynomialDerivative();
234
235 PolynomialFunction g0 = new PolynomialFunction(new double[] { k * (k + 1) });
236 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
237 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
238
239 PolynomialFunction Pk0g0 = Pk0.multiply(g0);
240 PolynomialFunction Pk1g1 = Pk1.multiply(g1);
241 PolynomialFunction Pk2g2 = Pk2.multiply(g2);
242
243 checkNullPolynomial(Pk0g0.add(Pk1g1.add(Pk2g2)));
244
245 }
246 }
247
248 @Test
249 void testLegendreOrthogonality() {
250 UnivariateFunction weight = new UnivariateFunction() {
251 @Override
252 public double value(double x) {
253 return 1;
254 }
255 };
256 for (int i = 0; i < 10; ++i) {
257 PolynomialFunction pi = PolynomialsUtils.createLegendrePolynomial(i);
258 for (int j = 0; j <= i; ++j) {
259 PolynomialFunction pj = PolynomialsUtils.createLegendrePolynomial(j);
260 checkOrthogonality(pi, pj, weight, -1, 1, 0.1, 1.0e-13);
261 }
262 }
263 }
264
265 @Test
266 void testHighDegreeLegendre() {
267 PolynomialsUtils.createLegendrePolynomial(40);
268 double[] l40 = PolynomialsUtils.createLegendrePolynomial(40).getCoefficients();
269 double denominator = 274877906944d;
270 double[] numerators = new double[] {
271 +34461632205d, -28258538408100d, +3847870979902950d, -207785032914759300d,
272 +5929294332103310025d, -103301483474866556880d, +1197358103913226000200d, -9763073770369381232400d,
273 +58171647881784229843050d, -260061484647976556945400d, +888315281771246239250340d, -2345767627188139419665400d,
274 +4819022625419112503443050d, -7710436200670580005508880d, +9566652323054238154983240d, -9104813935044723209570256d,
275 +6516550296251767619752905d, -3391858621221953912598660d, +1211378079007840683070950d, -265365894974690562152100d,
276 +26876802183334044115405d
277 };
278 for (int i = 0; i < l40.length; ++i) {
279 if (i % 2 == 0) {
280 double ci = numerators[i / 2] / denominator;
281 assertEquals(ci, l40[i], FastMath.abs(ci) * 1e-15);
282 } else {
283 assertEquals(0, l40[i], 0);
284 }
285 }
286 }
287
288 @Test
289 void testJacobiLegendre() {
290 for (int i = 0; i < 10; ++i) {
291 PolynomialFunction legendre = PolynomialsUtils.createLegendrePolynomial(i);
292 PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, 0, 0);
293 checkNullPolynomial(legendre.subtract(jacobi));
294 }
295 }
296
297 @Test
298 void testJacobiEvaluationAt1() {
299 for (int v = 0; v < 10; ++v) {
300 for (int w = 0; w < 10; ++w) {
301 for (int i = 0; i < 10; ++i) {
302 PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
303 double binomial = CombinatoricsUtils.binomialCoefficient(v + i, i);
304 assertTrue(Precision.equals(binomial, jacobi.value(1.0), 1));
305 }
306 }
307 }
308 }
309
310 @Test
311 void testJacobiOrthogonality() {
312 for (int v = 0; v < 5; ++v) {
313 for (int w = v; w < 5; ++w) {
314 final int vv = v;
315 final int ww = w;
316 UnivariateFunction weight = new UnivariateFunction() {
317 @Override
318 public double value(double x) {
319 return FastMath.pow(1 - x, vv) * FastMath.pow(1 + x, ww);
320 }
321 };
322 for (int i = 0; i < 10; ++i) {
323 PolynomialFunction pi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
324 for (int j = 0; j <= i; ++j) {
325 PolynomialFunction pj = PolynomialsUtils.createJacobiPolynomial(j, v, w);
326 checkOrthogonality(pi, pj, weight, -1, 1, 0.1, 1.0e-12);
327 }
328 }
329 }
330 }
331 }
332
333
334
335
336 @Test
337 void testJacobiKeyCoverage() {
338
339 final JacobiKey key = new JacobiKey(3, 4);
340
341 assertNotEquals(null, key);
342 assertNotEquals(key, new Object());
343 assertNotEquals(key, new JacobiKey(3, 5));
344 assertNotEquals(key, new JacobiKey(5, 4));
345 assertNotEquals(key, new JacobiKey(5, 5));
346 assertEquals(key, new JacobiKey(3, 4));
347 }
348
349 @Test
350 void testShift() {
351
352 PolynomialFunction f1x = new PolynomialFunction(new double[] { 1, 1, 2 });
353
354 PolynomialFunction f1x1
355 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 1));
356 checkPolynomial(f1x1, "4 + 5 x + 2 x^2");
357
358 PolynomialFunction f1xM1
359 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), -1));
360 checkPolynomial(f1xM1, "2 - 3 x + 2 x^2");
361
362 PolynomialFunction f1x3
363 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 3));
364 checkPolynomial(f1x3, "22 + 13 x + 2 x^2");
365
366
367 PolynomialFunction f2x = new PolynomialFunction(new double[]{2, 0, 3, 8, 0, 121});
368
369 PolynomialFunction f2x1
370 = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 1));
371 checkPolynomial(f2x1, "134 + 635 x + 1237 x^2 + 1218 x^3 + 605 x^4 + 121 x^5");
372
373 PolynomialFunction f2x3
374 = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 3));
375 checkPolynomial(f2x3, "29648 + 49239 x + 32745 x^2 + 10898 x^3 + 1815 x^4 + 121 x^5");
376 }
377
378
379 private void checkPolynomial(PolynomialFunction p, long denominator, String reference) {
380 PolynomialFunction q = new PolynomialFunction(new double[] { denominator});
381 assertEquals(reference, p.multiply(q).toString());
382 }
383
384 private void checkPolynomial(PolynomialFunction p, String reference) {
385 assertEquals(reference, p.toString());
386 }
387
388 private void checkNullPolynomial(PolynomialFunction p) {
389 for (double coefficient : p.getCoefficients()) {
390 assertEquals(0, coefficient, 1e-13);
391 }
392 }
393
394 private void checkOrthogonality(final PolynomialFunction p1,
395 final PolynomialFunction p2,
396 final UnivariateFunction weight,
397 final double a, final double b,
398 final double nonZeroThreshold,
399 final double zeroThreshold) {
400 UnivariateFunction f = new UnivariateFunction() {
401 @Override
402 public double value(double x) {
403 return weight.value(x) * p1.value(x) * p2.value(x);
404 }
405 };
406 double dotProduct =
407 new IterativeLegendreGaussIntegrator(5, 1.0e-9, 1.0e-8, 2, 15).integrate(1000000, f, a, b);
408 if (p1.degree() == p2.degree()) {
409
410 assertTrue(FastMath.abs(dotProduct) > nonZeroThreshold,
411 "I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct);
412 } else {
413
414 assertEquals(0.0, FastMath.abs(dotProduct), zeroThreshold, "I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct);
415 }
416 }
417 }