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