From 461c45aa7db5af1d78c369fff2f11293e5f31508 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Fri, 7 Oct 2011 16:31:04 +0000 Subject: [PATCH] Added Jacobi orthogonal polynomials. JIRA: MATH-687 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1180092 13f79535-47bb-0310-9956-ffa450edef68 --- .../polynomials/PolynomialsUtils.java | 106 +++++++++++++++++- src/site/xdoc/changes.xml | 3 + src/site/xdoc/userguide/analysis.xml | 2 +- .../polynomials/PolynomialsUtilsTest.java | 45 ++++++++ 4 files changed, 149 insertions(+), 7 deletions(-) diff --git a/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java b/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java index 6d8071bc2..e6e777587 100644 --- a/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java +++ b/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java @@ -17,6 +17,9 @@ package org.apache.commons.math.analysis.polynomials; import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; import org.apache.commons.math.fraction.BigFraction; import org.apache.commons.math.util.FastMath; @@ -31,16 +34,19 @@ import org.apache.commons.math.util.MathUtils; public class PolynomialsUtils { /** Coefficients for Chebyshev polynomials. */ - private static final ArrayList CHEBYSHEV_COEFFICIENTS; + private static final List CHEBYSHEV_COEFFICIENTS; /** Coefficients for Hermite polynomials. */ - private static final ArrayList HERMITE_COEFFICIENTS; + private static final List HERMITE_COEFFICIENTS; /** Coefficients for Laguerre polynomials. */ - private static final ArrayList LAGUERRE_COEFFICIENTS; + private static final List LAGUERRE_COEFFICIENTS; /** Coefficients for Legendre polynomials. */ - private static final ArrayList LEGENDRE_COEFFICIENTS; + private static final List LEGENDRE_COEFFICIENTS; + + /** Coefficients for Jacobi polynomials. */ + private static final Map> JACOBI_COEFFICIENTS; static { @@ -72,6 +78,9 @@ public class PolynomialsUtils { LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO); LEGENDRE_COEFFICIENTS.add(BigFraction.ONE); + // initialize map for Jacobi polynomials + JACOBI_COEFFICIENTS = new HashMap>(); + } /** @@ -185,6 +194,91 @@ public class PolynomialsUtils { }); } + /** + * Create a Jacobi polynomial. + *

Jacobi + * polynomials are orthogonal polynomials. + * They can be defined by the following recurrence relations: + * 

+     *        P0vw(X)   = 1
+     *        P-1vw(X)  = 0
+     *  2k(k + v + w)(2k + v + w - 2) Pkvw(X) = 
+     *  (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v2 - w2] Pk-1vw(X)
+     *  - 2(k + v - 1)(k + w - 1)(2k + v + w) Pk-2vw(X)
+     *

+ * @param degree degree of the polynomial + * @param v first exponent + * @param w second exponent + * @return Jacobi polynomial of specified degree + */ + public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) { + + // select the appropriate list + final JacobiKey key = new JacobiKey(v, w); + + if (!JACOBI_COEFFICIENTS.containsKey(key)) { + + // allocate a new list for v, w + final List list = new ArrayList(); + JACOBI_COEFFICIENTS.put(key, list); + + // Pv,w,0(x) = 1; + list.add(BigFraction.ONE); + + // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2 + list.add(new BigFraction(v - w, 2)); + list.add(new BigFraction(2 + v + w, 2)); + + } + + return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key), + new RecurrenceCoefficientsGenerator() { + /** {@inheritDoc} */ + public BigFraction[] generate(int k) { + k++; + final int kvw = k + v + w; + final int twoKvw = kvw + k; + final int twoKvwM1 = twoKvw - 1; + final int twoKvwM2 = twoKvw - 2; + final int den = 2 * k * kvw * twoKvwM2; + + return new BigFraction[] { + new BigFraction(twoKvwM1 * (v * v - w * w), den), + new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den), + new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den) + }; + } + }); + + } + + /** Inner class for Jacobi polynomials keys. */ + private static class JacobiKey { + + /** First exponent. */ + private final int v; + + /** Second exponent. */ + private final int w; + + /** Simple constructor. + * @param v first exponent + * @param w second exponent + */ + public JacobiKey(final int v, final int w) { + this.v = v; + this.w = w; + } + + /** Get hash code. + * @return hash code + */ + public int hashCode() { + return (v << 16) ^ w; + } + + } + /** * Compute the coefficients of the polynomial Ps(x) * whose values at point {@code x} will be the same as the those from the @@ -247,7 +341,7 @@ public class PolynomialsUtils { * @return coefficients array */ private static PolynomialFunction buildPolynomial(final int degree, - final ArrayList coefficients, + final List coefficients, final RecurrenceCoefficientsGenerator generator) { final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1; @@ -285,7 +379,7 @@ public class PolynomialsUtils { */ private static void computeUpToDegree(final int degree, final int maxDegree, final RecurrenceCoefficientsGenerator generator, - final ArrayList coefficients) { + final List coefficients) { int startK = (maxDegree - 1) * maxDegree / 2; for (int k = maxDegree; k < degree; ++k) { diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 22fb79695..757997e70 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,9 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added Jacobi polynomials. + Added "shift" method to compute the coefficients of a new polynomial whose values are the same as those of another polynomial but computed diff --git a/src/site/xdoc/userguide/analysis.xml b/src/site/xdoc/userguide/analysis.xml index c267b0c71..d43bf4ab3 100644 --- a/src/site/xdoc/userguide/analysis.xml +++ b/src/site/xdoc/userguide/analysis.xml @@ -522,7 +522,7 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY); coefficients arrays. The PolynomialsUtils utility class provides static factory methods to build - Chebyshev, Hermite, Lagrange and Legendre polynomials. Coefficients are + Chebyshev, Hermite, Jacobi, Laguerre and Legendre polynomials. Coefficients are computed using exact fractions so these factory methods can build polynomials up to any degree.

diff --git a/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java b/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java index a727e4636..336dd8d4d 100644 --- a/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java +++ b/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java @@ -19,6 +19,7 @@ package org.apache.commons.math.analysis.polynomials; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.analysis.integration.LegendreGaussIntegrator; import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; import org.junit.Assert; import org.junit.Test; @@ -272,6 +273,50 @@ public class PolynomialsUtilsTest { } } + @Test + public void testJacobiLegendre() { + for (int i = 0; i < 10; ++i) { + PolynomialFunction legendre = PolynomialsUtils.createLegendrePolynomial(i); + PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, 0, 0); + checkNullPolynomial(legendre.subtract(jacobi)); + } + } + + @Test + public void testJacobiEvaluationAt1() { + for (int v = 0; v < 10; ++v) { + for (int w = 0; w < 10; ++w) { + for (int i = 0; i < 10; ++i) { + PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w); + double binomial = MathUtils.binomialCoefficient(v + i, i); + Assert.assertTrue(MathUtils.equals(binomial, jacobi.value(1.0), 1)); + } + } + } + } + + @Test + public void testJacobiOrthogonality() { + for (int v = 0; v < 5; ++v) { + for (int w = v; w < 5; ++w) { + final int vv = v; + final int ww = w; + UnivariateRealFunction weight = new UnivariateRealFunction() { + public double value(double x) { + return FastMath.pow(1 - x, vv) * FastMath.pow(1 + x, ww); + } + }; + for (int i = 0; i < 10; ++i) { + PolynomialFunction pi = PolynomialsUtils.createJacobiPolynomial(i, v, w); + for (int j = 0; j <= i; ++j) { + PolynomialFunction pj = PolynomialsUtils.createJacobiPolynomial(j, v, w); + checkOrthogonality(pi, pj, weight, -1, 1, 0.1, 1.0e-12); + } + } + } + } + } + @Test public void testShift() { // f1(x) = 1 + x + 2 x^2