001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.polynomials;
018
019 import org.apache.commons.math3.analysis.UnivariateFunction;
020 import org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator;
021 import org.apache.commons.math3.util.ArithmeticUtils;
022 import org.apache.commons.math3.util.FastMath;
023 import org.apache.commons.math3.util.Precision;
024 import org.junit.Assert;
025 import org.junit.Test;
026
027 /**
028 * Tests the PolynomialsUtils class.
029 *
030 * @version $Id: PolynomialsUtilsTest.java 1377244 2012-08-25 10:05:13Z luc $
031 */
032 public class PolynomialsUtilsTest {
033
034 @Test
035 public void testFirstChebyshevPolynomials() {
036 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(3), "-3 x + 4 x^3");
037 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(2), "-1 + 2 x^2");
038 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(1), "x");
039 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(0), "1");
040
041 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(7), "-7 x + 56 x^3 - 112 x^5 + 64 x^7");
042 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(6), "-1 + 18 x^2 - 48 x^4 + 32 x^6");
043 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(5), "5 x - 20 x^3 + 16 x^5");
044 checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(4), "1 - 8 x^2 + 8 x^4");
045
046 }
047
048 @Test
049 public void testChebyshevBounds() {
050 for (int k = 0; k < 12; ++k) {
051 PolynomialFunction Tk = PolynomialsUtils.createChebyshevPolynomial(k);
052 for (double x = -1; x <= 1; x += 0.02) {
053 Assert.assertTrue(k + " " + Tk.value(x), FastMath.abs(Tk.value(x)) < (1 + 1e-12));
054 }
055 }
056 }
057
058 @Test
059 public void testChebyshevDifferentials() {
060 for (int k = 0; k < 12; ++k) {
061
062 PolynomialFunction Tk0 = PolynomialsUtils.createChebyshevPolynomial(k);
063 PolynomialFunction Tk1 = Tk0.polynomialDerivative();
064 PolynomialFunction Tk2 = Tk1.polynomialDerivative();
065
066 PolynomialFunction g0 = new PolynomialFunction(new double[] { k * k });
067 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -1});
068 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
069
070 PolynomialFunction Tk0g0 = Tk0.multiply(g0);
071 PolynomialFunction Tk1g1 = Tk1.multiply(g1);
072 PolynomialFunction Tk2g2 = Tk2.multiply(g2);
073
074 checkNullPolynomial(Tk0g0.add(Tk1g1.add(Tk2g2)));
075
076 }
077 }
078
079 @Test
080 public void testChebyshevOrthogonality() {
081 UnivariateFunction weight = new UnivariateFunction() {
082 public double value(double x) {
083 return 1 / FastMath.sqrt(1 - x * x);
084 }
085 };
086 for (int i = 0; i < 10; ++i) {
087 PolynomialFunction pi = PolynomialsUtils.createChebyshevPolynomial(i);
088 for (int j = 0; j <= i; ++j) {
089 PolynomialFunction pj = PolynomialsUtils.createChebyshevPolynomial(j);
090 checkOrthogonality(pi, pj, weight, -0.9999, 0.9999, 1.5, 0.03);
091 }
092 }
093 }
094
095 @Test
096 public void testFirstHermitePolynomials() {
097 checkPolynomial(PolynomialsUtils.createHermitePolynomial(3), "-12 x + 8 x^3");
098 checkPolynomial(PolynomialsUtils.createHermitePolynomial(2), "-2 + 4 x^2");
099 checkPolynomial(PolynomialsUtils.createHermitePolynomial(1), "2 x");
100 checkPolynomial(PolynomialsUtils.createHermitePolynomial(0), "1");
101
102 checkPolynomial(PolynomialsUtils.createHermitePolynomial(7), "-1680 x + 3360 x^3 - 1344 x^5 + 128 x^7");
103 checkPolynomial(PolynomialsUtils.createHermitePolynomial(6), "-120 + 720 x^2 - 480 x^4 + 64 x^6");
104 checkPolynomial(PolynomialsUtils.createHermitePolynomial(5), "120 x - 160 x^3 + 32 x^5");
105 checkPolynomial(PolynomialsUtils.createHermitePolynomial(4), "12 - 48 x^2 + 16 x^4");
106
107 }
108
109 @Test
110 public void testHermiteDifferentials() {
111 for (int k = 0; k < 12; ++k) {
112
113 PolynomialFunction Hk0 = PolynomialsUtils.createHermitePolynomial(k);
114 PolynomialFunction Hk1 = Hk0.polynomialDerivative();
115 PolynomialFunction Hk2 = Hk1.polynomialDerivative();
116
117 PolynomialFunction g0 = new PolynomialFunction(new double[] { 2 * k });
118 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
119 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1 });
120
121 PolynomialFunction Hk0g0 = Hk0.multiply(g0);
122 PolynomialFunction Hk1g1 = Hk1.multiply(g1);
123 PolynomialFunction Hk2g2 = Hk2.multiply(g2);
124
125 checkNullPolynomial(Hk0g0.add(Hk1g1.add(Hk2g2)));
126
127 }
128 }
129
130 @Test
131 public void testHermiteOrthogonality() {
132 UnivariateFunction weight = new UnivariateFunction() {
133 public double value(double x) {
134 return FastMath.exp(-x * x);
135 }
136 };
137 for (int i = 0; i < 10; ++i) {
138 PolynomialFunction pi = PolynomialsUtils.createHermitePolynomial(i);
139 for (int j = 0; j <= i; ++j) {
140 PolynomialFunction pj = PolynomialsUtils.createHermitePolynomial(j);
141 checkOrthogonality(pi, pj, weight, -50, 50, 1.5, 1.0e-8);
142 }
143 }
144 }
145
146 @Test
147 public void testFirstLaguerrePolynomials() {
148 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(3), 6l, "6 - 18 x + 9 x^2 - x^3");
149 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(2), 2l, "2 - 4 x + x^2");
150 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(1), 1l, "1 - x");
151 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(0), 1l, "1");
152
153 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(7), 5040l,
154 "5040 - 35280 x + 52920 x^2 - 29400 x^3"
155 + " + 7350 x^4 - 882 x^5 + 49 x^6 - x^7");
156 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(6), 720l,
157 "720 - 4320 x + 5400 x^2 - 2400 x^3 + 450 x^4"
158 + " - 36 x^5 + x^6");
159 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(5), 120l,
160 "120 - 600 x + 600 x^2 - 200 x^3 + 25 x^4 - x^5");
161 checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(4), 24l,
162 "24 - 96 x + 72 x^2 - 16 x^3 + x^4");
163
164 }
165
166 @Test
167 public void testLaguerreDifferentials() {
168 for (int k = 0; k < 12; ++k) {
169
170 PolynomialFunction Lk0 = PolynomialsUtils.createLaguerrePolynomial(k);
171 PolynomialFunction Lk1 = Lk0.polynomialDerivative();
172 PolynomialFunction Lk2 = Lk1.polynomialDerivative();
173
174 PolynomialFunction g0 = new PolynomialFunction(new double[] { k });
175 PolynomialFunction g1 = new PolynomialFunction(new double[] { 1, -1 });
176 PolynomialFunction g2 = new PolynomialFunction(new double[] { 0, 1 });
177
178 PolynomialFunction Lk0g0 = Lk0.multiply(g0);
179 PolynomialFunction Lk1g1 = Lk1.multiply(g1);
180 PolynomialFunction Lk2g2 = Lk2.multiply(g2);
181
182 checkNullPolynomial(Lk0g0.add(Lk1g1.add(Lk2g2)));
183
184 }
185 }
186
187 @Test
188 public void testLaguerreOrthogonality() {
189 UnivariateFunction weight = new UnivariateFunction() {
190 public double value(double x) {
191 return FastMath.exp(-x);
192 }
193 };
194 for (int i = 0; i < 10; ++i) {
195 PolynomialFunction pi = PolynomialsUtils.createLaguerrePolynomial(i);
196 for (int j = 0; j <= i; ++j) {
197 PolynomialFunction pj = PolynomialsUtils.createLaguerrePolynomial(j);
198 checkOrthogonality(pi, pj, weight, 0.0, 100.0, 0.99999, 1.0e-13);
199 }
200 }
201 }
202
203 @Test
204 public void testFirstLegendrePolynomials() {
205 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(3), 2l, "-3 x + 5 x^3");
206 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(2), 2l, "-1 + 3 x^2");
207 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(1), 1l, "x");
208 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(0), 1l, "1");
209
210 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(7), 16l, "-35 x + 315 x^3 - 693 x^5 + 429 x^7");
211 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(6), 16l, "-5 + 105 x^2 - 315 x^4 + 231 x^6");
212 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(5), 8l, "15 x - 70 x^3 + 63 x^5");
213 checkPolynomial(PolynomialsUtils.createLegendrePolynomial(4), 8l, "3 - 30 x^2 + 35 x^4");
214
215 }
216
217 @Test
218 public void testLegendreDifferentials() {
219 for (int k = 0; k < 12; ++k) {
220
221 PolynomialFunction Pk0 = PolynomialsUtils.createLegendrePolynomial(k);
222 PolynomialFunction Pk1 = Pk0.polynomialDerivative();
223 PolynomialFunction Pk2 = Pk1.polynomialDerivative();
224
225 PolynomialFunction g0 = new PolynomialFunction(new double[] { k * (k + 1) });
226 PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
227 PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
228
229 PolynomialFunction Pk0g0 = Pk0.multiply(g0);
230 PolynomialFunction Pk1g1 = Pk1.multiply(g1);
231 PolynomialFunction Pk2g2 = Pk2.multiply(g2);
232
233 checkNullPolynomial(Pk0g0.add(Pk1g1.add(Pk2g2)));
234
235 }
236 }
237
238 @Test
239 public void testLegendreOrthogonality() {
240 UnivariateFunction weight = new UnivariateFunction() {
241 public double value(double x) {
242 return 1;
243 }
244 };
245 for (int i = 0; i < 10; ++i) {
246 PolynomialFunction pi = PolynomialsUtils.createLegendrePolynomial(i);
247 for (int j = 0; j <= i; ++j) {
248 PolynomialFunction pj = PolynomialsUtils.createLegendrePolynomial(j);
249 checkOrthogonality(pi, pj, weight, -1, 1, 0.1, 1.0e-13);
250 }
251 }
252 }
253
254 @Test
255 public void testHighDegreeLegendre() {
256 PolynomialsUtils.createLegendrePolynomial(40);
257 double[] l40 = PolynomialsUtils.createLegendrePolynomial(40).getCoefficients();
258 double denominator = 274877906944d;
259 double[] numerators = new double[] {
260 +34461632205d, -28258538408100d, +3847870979902950d, -207785032914759300d,
261 +5929294332103310025d, -103301483474866556880d, +1197358103913226000200d, -9763073770369381232400d,
262 +58171647881784229843050d, -260061484647976556945400d, +888315281771246239250340d, -2345767627188139419665400d,
263 +4819022625419112503443050d, -7710436200670580005508880d, +9566652323054238154983240d, -9104813935044723209570256d,
264 +6516550296251767619752905d, -3391858621221953912598660d, +1211378079007840683070950d, -265365894974690562152100d,
265 +26876802183334044115405d
266 };
267 for (int i = 0; i < l40.length; ++i) {
268 if (i % 2 == 0) {
269 double ci = numerators[i / 2] / denominator;
270 Assert.assertEquals(ci, l40[i], FastMath.abs(ci) * 1e-15);
271 } else {
272 Assert.assertEquals(0, l40[i], 0);
273 }
274 }
275 }
276
277 @Test
278 public void testJacobiLegendre() {
279 for (int i = 0; i < 10; ++i) {
280 PolynomialFunction legendre = PolynomialsUtils.createLegendrePolynomial(i);
281 PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, 0, 0);
282 checkNullPolynomial(legendre.subtract(jacobi));
283 }
284 }
285
286 @Test
287 public void testJacobiEvaluationAt1() {
288 for (int v = 0; v < 10; ++v) {
289 for (int w = 0; w < 10; ++w) {
290 for (int i = 0; i < 10; ++i) {
291 PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
292 double binomial = ArithmeticUtils.binomialCoefficient(v + i, i);
293 Assert.assertTrue(Precision.equals(binomial, jacobi.value(1.0), 1));
294 }
295 }
296 }
297 }
298
299 @Test
300 public void testJacobiOrthogonality() {
301 for (int v = 0; v < 5; ++v) {
302 for (int w = v; w < 5; ++w) {
303 final int vv = v;
304 final int ww = w;
305 UnivariateFunction weight = new UnivariateFunction() {
306 public double value(double x) {
307 return FastMath.pow(1 - x, vv) * FastMath.pow(1 + x, ww);
308 }
309 };
310 for (int i = 0; i < 10; ++i) {
311 PolynomialFunction pi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
312 for (int j = 0; j <= i; ++j) {
313 PolynomialFunction pj = PolynomialsUtils.createJacobiPolynomial(j, v, w);
314 checkOrthogonality(pi, pj, weight, -1, 1, 0.1, 1.0e-12);
315 }
316 }
317 }
318 }
319 }
320
321 @Test
322 public void testShift() {
323 // f1(x) = 1 + x + 2 x^2
324 PolynomialFunction f1x = new PolynomialFunction(new double[] { 1, 1, 2 });
325
326 PolynomialFunction f1x1
327 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 1));
328 checkPolynomial(f1x1, "4 + 5 x + 2 x^2");
329
330 PolynomialFunction f1xM1
331 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), -1));
332 checkPolynomial(f1xM1, "2 - 3 x + 2 x^2");
333
334 PolynomialFunction f1x3
335 = new PolynomialFunction(PolynomialsUtils.shift(f1x.getCoefficients(), 3));
336 checkPolynomial(f1x3, "22 + 13 x + 2 x^2");
337
338 // f2(x) = 2 + 3 x^2 + 8 x^3 + 121 x^5
339 PolynomialFunction f2x = new PolynomialFunction(new double[]{2, 0, 3, 8, 0, 121});
340
341 PolynomialFunction f2x1
342 = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 1));
343 checkPolynomial(f2x1, "134 + 635 x + 1237 x^2 + 1218 x^3 + 605 x^4 + 121 x^5");
344
345 PolynomialFunction f2x3
346 = new PolynomialFunction(PolynomialsUtils.shift(f2x.getCoefficients(), 3));
347 checkPolynomial(f2x3, "29648 + 49239 x + 32745 x^2 + 10898 x^3 + 1815 x^4 + 121 x^5");
348 }
349
350
351 private void checkPolynomial(PolynomialFunction p, long denominator, String reference) {
352 PolynomialFunction q = new PolynomialFunction(new double[] { denominator});
353 Assert.assertEquals(reference, p.multiply(q).toString());
354 }
355
356 private void checkPolynomial(PolynomialFunction p, String reference) {
357 Assert.assertEquals(reference, p.toString());
358 }
359
360 private void checkNullPolynomial(PolynomialFunction p) {
361 for (double coefficient : p.getCoefficients()) {
362 Assert.assertEquals(0, coefficient, 1e-13);
363 }
364 }
365
366 private void checkOrthogonality(final PolynomialFunction p1,
367 final PolynomialFunction p2,
368 final UnivariateFunction weight,
369 final double a, final double b,
370 final double nonZeroThreshold,
371 final double zeroThreshold) {
372 UnivariateFunction f = new UnivariateFunction() {
373 public double value(double x) {
374 return weight.value(x) * p1.value(x) * p2.value(x);
375 }
376 };
377 double dotProduct =
378 new IterativeLegendreGaussIntegrator(5, 1.0e-9, 1.0e-8, 2, 15).integrate(1000000, f, a, b);
379 if (p1.degree() == p2.degree()) {
380 // integral should be non-zero
381 Assert.assertTrue("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
382 FastMath.abs(dotProduct) > nonZeroThreshold);
383 } else {
384 // integral should be zero
385 Assert.assertEquals("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
386 0.0, FastMath.abs(dotProduct), zeroThreshold);
387 }
388 }
389 }