commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1180092 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/analysis/polynomials/ site/xdoc/ site/xdoc/userguide/ test/java/org/apache/commons/math/analysis/polynomials/
Date Fri, 07 Oct 2011 16:31:04 GMT
Author: luc
Date: Fri Oct  7 16:31:04 2011
New Revision: 1180092

URL: http://svn.apache.org/viewvc?rev=1180092&view=rev
Log:
Added Jacobi orthogonal polynomials.

JIRA: MATH-687

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java?rev=1180092&r1=1180091&r2=1180092&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
Fri Oct  7 16:31:04 2011
@@ -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.Math
 public class PolynomialsUtils {
 
     /** Coefficients for Chebyshev polynomials. */
-    private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
+    private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
 
     /** Coefficients for Hermite polynomials. */
-    private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
+    private static final List<BigFraction> HERMITE_COEFFICIENTS;
 
     /** Coefficients for Laguerre polynomials. */
-    private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
+    private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
 
     /** Coefficients for Legendre polynomials. */
-    private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
+    private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
+
+    /** Coefficients for Jacobi polynomials. */
+    private static final Map<JacobiKey, List<BigFraction>> 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<JacobiKey, List<BigFraction>>();
+
     }
 
     /**
@@ -186,6 +195,91 @@ public class PolynomialsUtils {
     }
 
     /**
+     * Create a Jacobi polynomial.
+     * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi

+     * polynomials</a> are orthogonal polynomials.
+     * They can be defined by the following recurrence relations:
+     * <pre>
+     *        P<sub>0</sub><sup>vw</sup>(X)   = 1
+     *        P<sub>-1</sub><sup>vw</sup>(X)  = 0
+     *  2k(k + v + w)(2k + v + w - 2) P<sub>k</sub><sup>vw</sup>(X)
= 
+     *  (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v<sup>2</sup> - w<sup>2</sup>]
P<sub>k-1</sub><sup>vw</sup>(X)
+     *  - 2(k + v - 1)(k + w - 1)(2k + v + w) P<sub>k-2</sub><sup>vw</sup>(X)
+     * </pre></p>
+     * @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<BigFraction> list = new ArrayList<BigFraction>();
+            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 <code>P<sub>s</sub>(x)</code>
      * whose values at point {@code x} will be the same as the those from the
      * original polynomial <code>P(x)</code> when computed at {@code x + shift}.
@@ -247,7 +341,7 @@ public class PolynomialsUtils {
      * @return coefficients array
      */
     private static PolynomialFunction buildPolynomial(final int degree,
-                                                      final ArrayList<BigFraction>
coefficients,
+                                                      final List<BigFraction> 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<BigFraction> coefficients)
{
+                                          final List<BigFraction> coefficients) {
 
         int startK = (maxDegree - 1) * maxDegree / 2;
         for (int k = maxDegree; k < degree; ++k) {

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1180092&r1=1180091&r2=1180092&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Oct  7 16:31:04 2011
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-687" due-to="Romain di Costanzo">
+         Added Jacobi polynomials.
+      </action>
       <action dev="erans" type="add" issue="MATH-683" due-to="Romain di Costanzo">
          Added "shift" method to compute the coefficients of a new polynomial
          whose values are the same as those of another polynomial but computed

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml?rev=1180092&r1=1180091&r2=1180092&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml Fri Oct  7 16:31:04 2011
@@ -522,7 +522,7 @@ System.out println("f(" + interpolationX
           coefficients arrays. The
           <a href="../apidocs/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.html">
           PolynomialsUtils</a> 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.
         </p>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java?rev=1180092&r1=1180091&r2=1180092&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
Fri Oct  7 16:31:04 2011
@@ -19,6 +19,7 @@ package org.apache.commons.math.analysis
 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;
 
@@ -273,6 +274,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
         PolynomialFunction f1x = new PolynomialFunction(new double[] { 1, 1, 2 });



Mime
View raw message