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