View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
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   * Tests the PolynomialsUtils class.
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     @Test
331     public void testShift() {
332         // f1(x) = 1 + x + 2 x^2
333         PolynomialFunction f1x = new PolynomialFunction(new double[] { 1, 1, 2 });
334 
335         PolynomialFunction f1x1
336             = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 1));
337         checkPolynomial(f1x1, "4 + 5 x + 2 x^2");
338 
339         PolynomialFunction f1xM1
340             = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), -1));
341         checkPolynomial(f1xM1, "2 - 3 x + 2 x^2");
342 
343         PolynomialFunction f1x3
344             = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 3));
345         checkPolynomial(f1x3, "22 + 13 x + 2 x^2");
346 
347         // f2(x) = 2 + 3 x^2 + 8 x^3 + 121 x^5
348         PolynomialFunction f2x = new PolynomialFunction(new double[]{2, 0, 3, 8, 0, 121});
349 
350         PolynomialFunction f2x1
351             = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 1));
352         checkPolynomial(f2x1, "134 + 635 x + 1237 x^2 + 1218 x^3 + 605 x^4 + 121 x^5");
353 
354         PolynomialFunction f2x3
355             = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 3));
356         checkPolynomial(f2x3, "29648 + 49239 x + 32745 x^2 + 10898 x^3 + 1815 x^4 + 121 x^5");
357     }
358 
359 
360     private void checkPolynomial(PolynomialFunction p, long denominator, String reference) {
361         PolynomialFunction q = new PolynomialFunction(new double[] { denominator});
362         Assert.assertEquals(reference, p.multiply(q).toString());
363     }
364 
365     private void checkPolynomial(PolynomialFunction p, String reference) {
366         Assert.assertEquals(reference, p.toString());
367     }
368 
369     private void checkNullPolynomial(PolynomialFunction p) {
370         for (double coefficient : p.getCoefficients()) {
371             Assert.assertEquals(0, coefficient, 1e-13);
372         }
373     }
374 
375     private void checkOrthogonality(final PolynomialFunction p1,
376                                     final PolynomialFunction p2,
377                                     final UnivariateFunction weight,
378                                     final double a, final double b,
379                                     final double nonZeroThreshold,
380                                     final double zeroThreshold) {
381         UnivariateFunction f = new UnivariateFunction() {
382             @Override
383             public double value(double x) {
384                 return weight.value(x) * p1.value(x) * p2.value(x);
385             }
386         };
387         double dotProduct =
388                 new IterativeLegendreGaussIntegrator(5, 1.0e-9, 1.0e-8, 2, 15).integrate(1000000, f, a, b);
389         if (p1.degree() == p2.degree()) {
390             // integral should be non-zero
391             Assert.assertTrue("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
392                               FastMath.abs(dotProduct) > nonZeroThreshold);
393         } else {
394             // integral should be zero
395             Assert.assertEquals("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
396                                 0.0, FastMath.abs(dotProduct), zeroThreshold);
397         }
398     }
399 }