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
This commit is contained in:
Luc Maisonobe
parent
ca133d8940
commit
461c45aa7d
|
|
@ -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<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>>();
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -185,6 +194,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
|
||||
|
|
@ -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) {
|
||||
|
|
|
|||
|
|
@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
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
|
||||
|
|
|
|||
|
|
@ -522,7 +522,7 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
|
|||
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>
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Reference in New Issue