You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by lu...@apache.org on 2007/02/12 20:27:17 UTC
svn commit: r506592 - in /jakarta/commons/proper/math/trunk/src:
java/org/apache/commons/math/ java/org/apache/commons/math/analysis/
test/org/apache/commons/math/
Author: luc
Date: Mon Feb 12 11:27:16 2007
New Revision: 506592
URL: http://svn.apache.org/viewvc?view=rev&rev=506592
Log:
Added and used a specialized exception for duplicate abscissas in sampled functions
Added:
jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java (with props)
jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java (with props)
Modified:
jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java (contents, props changed)
jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java (contents, props changed)
Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java?view=auto&rev=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java Mon Feb 12 11:27:16 2007
@@ -0,0 +1,47 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math;
+
+/**
+ * Exeption thrown when a sample contains several entries at the same abscissa.
+ * @version $Revision:$
+ */
+public class DuplicateSampleAbscissaException extends MathException {
+
+ /** Serializable version identifier */
+ private static final long serialVersionUID = -2271007547170169872L;
+
+ /**
+ * Construct an exception indicating the duplicate abscissa.
+ * @param abscissa duplicate abscissa
+ * @param i1 index of one entry having the duplicate abscissa
+ * @param i2 index of another entry having the duplicate abscissa
+ */
+ public DuplicateSampleAbscissaException(double abscissa, int i1, int i2) {
+ super("Abscissa {0} is duplicated at both indices {1} and {2}",
+ new Object[] { new Double(abscissa), new Integer(i1), new Integer(i2) });
+ }
+
+ /**
+ * Get the duplicate abscissa.
+ * @return duplicate abscissa
+ */
+ public double getDuplicateAbscissa() {
+ return ((Double) getArguments()[0]).doubleValue();
+ }
+
+}
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
------------------------------------------------------------------------------
svn:eol-style = native
Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java?view=diff&rev=506592&r1=506591&r2=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java (original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java Mon Feb 12 11:27:16 2007
@@ -1,122 +1,122 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.math.analysis;
-
-import java.io.Serializable;
-import org.apache.commons.math.MathException;
-
-/**
- * Implements the <a href="
- * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
- * Divided Difference Algorithm</a> for interpolation of real univariate
- * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
- * ISBN 038795452X, chapter 2.
- * <p>
- * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
- * this class provides an easy-to-use interface to it.
- *
- * @version $Revision$ $Date$
- */
-public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
- Serializable {
-
- /** serializable version identifier */
- static final long serialVersionUID = 107049519551235069L;
-
- /**
- * Computes an interpolating function for the data set.
- *
- * @param x the interpolating points array
- * @param y the interpolating values array
- * @return a function which interpolates the data set
- * @throws MathException if arguments are invalid
- */
- public UnivariateRealFunction interpolate(double x[], double y[]) throws
- MathException {
-
- /**
- * a[] and c[] are defined in the general formula of Newton form:
- * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
- * a[n](x-c[0])(x-c[1])...(x-c[n-1])
- */
- double a[], c[];
-
- PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
-
- /**
- * When used for interpolation, the Newton form formula becomes
- * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
- * f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
- * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
- * <p>
- * Note x[], y[], a[] have the same length but c[]'s size is one less.
- */
- c = new double[x.length-1];
- for (int i = 0; i < c.length; i++) {
- c[i] = x[i];
- }
- a = computeDividedDifference(x, y);
-
- PolynomialFunctionNewtonForm p;
- p = new PolynomialFunctionNewtonForm(a, c);
- return p;
- }
-
- /**
- * Returns a copy of the divided difference array.
- * <p>
- * The divided difference array is defined recursively by <pre>
- * f[x0] = f(x0)
- * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
- * </pre><p>
- * The computational complexity is O(N^2).
- *
- * @return a fresh copy of the divided difference array
- * @throws MathException if any abscissas coincide
- */
- protected static double[] computeDividedDifference(double x[], double y[])
- throws MathException {
-
- int i, j, n;
- double divdiff[], a[], denominator;
-
- PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
-
- n = x.length;
- divdiff = new double[n];
- for (i = 0; i < n; i++) {
- divdiff[i] = y[i]; // initialization
- }
-
- a = new double [n];
- a[0] = divdiff[0];
- for (i = 1; i < n; i++) {
- for (j = 0; j < n-i; j++) {
- denominator = x[j+i] - x[j];
- if (denominator == 0.0) {
- // This happens only when two abscissas are identical.
- throw new MathException
- ("Identical abscissas cause division by zero.");
- }
- divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
- }
- a[i] = divdiff[0];
- }
-
- return a;
- }
-}
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.DuplicateSampleAbscissaException;
+
+/**
+ * Implements the <a href="
+ * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
+ * Divided Difference Algorithm</a> for interpolation of real univariate
+ * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
+ * ISBN 038795452X, chapter 2.
+ * <p>
+ * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
+ * this class provides an easy-to-use interface to it.
+ *
+ * @version $Revision$ $Date$
+ */
+public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
+ Serializable {
+
+ /** serializable version identifier */
+ private static final long serialVersionUID = 107049519551235069L;
+
+ /**
+ * Computes an interpolating function for the data set.
+ *
+ * @param x the interpolating points array
+ * @param y the interpolating values array
+ * @return a function which interpolates the data set
+ * @throws DuplicateSampleAbscissaException if arguments are invalid
+ */
+ public UnivariateRealFunction interpolate(double x[], double y[]) throws
+ DuplicateSampleAbscissaException {
+
+ /**
+ * a[] and c[] are defined in the general formula of Newton form:
+ * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
+ * a[n](x-c[0])(x-c[1])...(x-c[n-1])
+ */
+ double a[], c[];
+
+ PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
+
+ /**
+ * When used for interpolation, the Newton form formula becomes
+ * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
+ * f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
+ * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
+ * <p>
+ * Note x[], y[], a[] have the same length but c[]'s size is one less.
+ */
+ c = new double[x.length-1];
+ for (int i = 0; i < c.length; i++) {
+ c[i] = x[i];
+ }
+ a = computeDividedDifference(x, y);
+
+ PolynomialFunctionNewtonForm p;
+ p = new PolynomialFunctionNewtonForm(a, c);
+ return p;
+ }
+
+ /**
+ * Returns a copy of the divided difference array.
+ * <p>
+ * The divided difference array is defined recursively by <pre>
+ * f[x0] = f(x0)
+ * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
+ * </pre><p>
+ * The computational complexity is O(N^2).
+ *
+ * @return a fresh copy of the divided difference array
+ * @throws DuplicateSampleAbscissaException if any abscissas coincide
+ */
+ protected static double[] computeDividedDifference(double x[], double y[])
+ throws DuplicateSampleAbscissaException {
+
+ int i, j, n;
+ double divdiff[], a[], denominator;
+
+ PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
+
+ n = x.length;
+ divdiff = new double[n];
+ for (i = 0; i < n; i++) {
+ divdiff[i] = y[i]; // initialization
+ }
+
+ a = new double [n];
+ a[0] = divdiff[0];
+ for (i = 1; i < n; i++) {
+ for (j = 0; j < n-i; j++) {
+ denominator = x[j+i] - x[j];
+ if (denominator == 0.0) {
+ // This happens only when two abscissas are identical.
+ throw new DuplicateSampleAbscissaException(x[j], j, j+i);
+ }
+ divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
+ }
+ a[i] = divdiff[0];
+ }
+
+ return a;
+ }
+}
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
('svn:executable' removed)
Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java?view=diff&rev=506592&r1=506591&r2=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java (original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java Mon Feb 12 11:27:16 2007
@@ -1,291 +1,295 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.math.analysis;
-
-import java.io.Serializable;
-import org.apache.commons.math.FunctionEvaluationException;
-
-/**
- * Implements the representation of a real polynomial function in
- * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
- * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
- * Analysis</b>, ISBN 038795452X, chapter 2.
- * <p>
- * The approximated function should be smooth enough for Lagrange polynomial
- * to work well. Otherwise, consider using splines instead.
- *
- * @version $Revision$ $Date$
- */
-public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
- Serializable {
-
- /** serializable version identifier */
- static final long serialVersionUID = -3965199246151093920L;
-
- /**
- * The coefficients of the polynomial, ordered by degree -- i.e.
- * coefficients[0] is the constant term and coefficients[n] is the
- * coefficient of x^n where n is the degree of the polynomial.
- */
- private double coefficients[];
-
- /**
- * Interpolating points (abscissas) and the function values at these points.
- */
- private double x[], y[];
-
- /**
- * Whether the polynomial coefficients are available.
- */
- private boolean coefficientsComputed;
-
- /**
- * Construct a Lagrange polynomial with the given abscissas and function
- * values. The order of interpolating points are not important.
- * <p>
- * The constructor makes copy of the input arrays and assigns them.
- *
- * @param x interpolating points
- * @param y function values at interpolating points
- * @throws IllegalArgumentException if input arrays are not valid
- */
- PolynomialFunctionLagrangeForm(double x[], double y[]) throws
- IllegalArgumentException {
-
- verifyInterpolationArray(x, y);
- this.x = new double[x.length];
- this.y = new double[y.length];
- System.arraycopy(x, 0, this.x, 0, x.length);
- System.arraycopy(y, 0, this.y, 0, y.length);
- coefficientsComputed = false;
- }
-
- /**
- * Calculate the function value at the given point.
- *
- * @param z the point at which the function value is to be computed
- * @return the function value
- * @throws FunctionEvaluationException if a runtime error occurs
- * @see UnivariateRealFunction#value(double)
- */
- public double value(double z) throws FunctionEvaluationException {
- return evaluate(x, y, z);
- }
-
- /**
- * Returns the degree of the polynomial.
- *
- * @return the degree of the polynomial
- */
- public int degree() {
- return x.length - 1;
- }
-
- /**
- * Returns a copy of the interpolating points array.
- * <p>
- * Changes made to the returned copy will not affect the polynomial.
- *
- * @return a fresh copy of the interpolating points array
- */
- public double[] getInterpolatingPoints() {
- double[] out = new double[x.length];
- System.arraycopy(x, 0, out, 0, x.length);
- return out;
- }
-
- /**
- * Returns a copy of the interpolating values array.
- * <p>
- * Changes made to the returned copy will not affect the polynomial.
- *
- * @return a fresh copy of the interpolating values array
- */
- public double[] getInterpolatingValues() {
- double[] out = new double[y.length];
- System.arraycopy(y, 0, out, 0, y.length);
- return out;
- }
-
- /**
- * Returns a copy of the coefficients array.
- * <p>
- * Changes made to the returned copy will not affect the polynomial.
- *
- * @return a fresh copy of the coefficients array
- */
- public double[] getCoefficients() {
- if (!coefficientsComputed) {
- computeCoefficients();
- }
- double[] out = new double[coefficients.length];
- System.arraycopy(coefficients, 0, out, 0, coefficients.length);
- return out;
- }
-
- /**
- * Evaluate the Lagrange polynomial using
- * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
- * Neville's Algorithm</a>. It takes O(N^2) time.
- * <p>
- * This function is made public static so that users can call it directly
- * without instantiating PolynomialFunctionLagrangeForm object.
- *
- * @param x the interpolating points array
- * @param y the interpolating values array
- * @param z the point at which the function value is to be computed
- * @return the function value
- * @throws FunctionEvaluationException if a runtime error occurs
- * @throws IllegalArgumentException if inputs are not valid
- */
- public static double evaluate(double x[], double y[], double z) throws
- FunctionEvaluationException, IllegalArgumentException {
-
- int i, j, n, nearest = 0;
- double value, c[], d[], tc, td, divider, w, dist, min_dist;
-
- verifyInterpolationArray(x, y);
-
- n = x.length;
- c = new double[n];
- d = new double[n];
- min_dist = Double.POSITIVE_INFINITY;
- for (i = 0; i < n; i++) {
- // initialize the difference arrays
- c[i] = y[i];
- d[i] = y[i];
- // find out the abscissa closest to z
- dist = Math.abs(z - x[i]);
- if (dist < min_dist) {
- nearest = i;
- min_dist = dist;
- }
- }
-
- // initial approximation to the function value at z
- value = y[nearest];
-
- for (i = 1; i < n; i++) {
- for (j = 0; j < n-i; j++) {
- tc = x[j] - z;
- td = x[i+j] - z;
- divider = x[j] - x[i+j];
- if (divider == 0.0) {
- // This happens only when two abscissas are identical.
- throw new FunctionEvaluationException(z,
- "Identical abscissas cause division by zero: x[" +
- i + "] = x[" + (i+j) + "] = " + x[i]);
- }
- // update the difference arrays
- w = (c[j+1] - d[j]) / divider;
- c[j] = tc * w;
- d[j] = td * w;
- }
- // sum up the difference terms to get the final value
- if (nearest < 0.5*(n-i+1)) {
- value += c[nearest]; // fork down
- } else {
- nearest--;
- value += d[nearest]; // fork up
- }
- }
-
- return value;
- }
-
- /**
- * Calculate the coefficients of Lagrange polynomial from the
- * interpolation data. It takes O(N^2) time.
- * <p>
- * Note this computation can be ill-conditioned. Use with caution
- * and only when it is necessary.
- *
- * @throws ArithmeticException if any abscissas coincide
- */
- protected void computeCoefficients() throws ArithmeticException {
- int i, j, n;
- double c[], tc[], d, t;
-
- n = degree() + 1;
- coefficients = new double[n];
- for (i = 0; i < n; i++) {
- coefficients[i] = 0.0;
- }
-
- // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
- c = new double[n+1];
- c[0] = 1.0;
- for (i = 0; i < n; i++) {
- for (j = i; j > 0; j--) {
- c[j] = c[j-1] - c[j] * x[i];
- }
- c[0] *= (-x[i]);
- c[i+1] = 1;
- }
-
- tc = new double[n];
- for (i = 0; i < n; i++) {
- // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
- d = 1;
- for (j = 0; j < n; j++) {
- if (i != j) {
- d *= (x[i] - x[j]);
- }
- }
- if (d == 0.0) {
- // This happens only when two abscissas are identical.
- throw new ArithmeticException
- ("Identical abscissas cause division by zero.");
- }
- t = y[i] / d;
- // Lagrange polynomial is the sum of n terms, each of which is a
- // polynomial of degree n-1. tc[] are the coefficients of the i-th
- // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
- tc[n-1] = c[n]; // actually c[n] = 1
- coefficients[n-1] += t * tc[n-1];
- for (j = n-2; j >= 0; j--) {
- tc[j] = c[j+1] + tc[j+1] * x[i];
- coefficients[j] += t * tc[j];
- }
- }
-
- coefficientsComputed = true;
- }
-
- /**
- * Verifies that the interpolation arrays are valid.
- * <p>
- * The interpolating points must be distinct. However it is not
- * verified here, it is checked in evaluate() and computeCoefficients().
- *
- * @throws IllegalArgumentException if not valid
- * @see #evaluate(double[], double[], double)
- * @see #computeCoefficients()
- */
- protected static void verifyInterpolationArray(double x[], double y[]) throws
- IllegalArgumentException {
-
- if (x.length < 2 || y.length < 2) {
- throw new IllegalArgumentException
- ("Interpolation requires at least two points.");
- }
- if (x.length != y.length) {
- throw new IllegalArgumentException
- ("Abscissa and value arrays must have the same length.");
- }
- }
-}
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.DuplicateSampleAbscissaException;
+import org.apache.commons.math.FunctionEvaluationException;
+
+/**
+ * Implements the representation of a real polynomial function in
+ * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
+ * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
+ * Analysis</b>, ISBN 038795452X, chapter 2.
+ * <p>
+ * The approximated function should be smooth enough for Lagrange polynomial
+ * to work well. Otherwise, consider using splines instead.
+ *
+ * @version $Revision$ $Date$
+ */
+public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
+ Serializable {
+
+ /** serializable version identifier */
+ static final long serialVersionUID = -3965199246151093920L;
+
+ /**
+ * The coefficients of the polynomial, ordered by degree -- i.e.
+ * coefficients[0] is the constant term and coefficients[n] is the
+ * coefficient of x^n where n is the degree of the polynomial.
+ */
+ private double coefficients[];
+
+ /**
+ * Interpolating points (abscissas) and the function values at these points.
+ */
+ private double x[], y[];
+
+ /**
+ * Whether the polynomial coefficients are available.
+ */
+ private boolean coefficientsComputed;
+
+ /**
+ * Construct a Lagrange polynomial with the given abscissas and function
+ * values. The order of interpolating points are not important.
+ * <p>
+ * The constructor makes copy of the input arrays and assigns them.
+ *
+ * @param x interpolating points
+ * @param y function values at interpolating points
+ * @throws IllegalArgumentException if input arrays are not valid
+ */
+ PolynomialFunctionLagrangeForm(double x[], double y[]) throws
+ IllegalArgumentException {
+
+ verifyInterpolationArray(x, y);
+ this.x = new double[x.length];
+ this.y = new double[y.length];
+ System.arraycopy(x, 0, this.x, 0, x.length);
+ System.arraycopy(y, 0, this.y, 0, y.length);
+ coefficientsComputed = false;
+ }
+
+ /**
+ * Calculate the function value at the given point.
+ *
+ * @param z the point at which the function value is to be computed
+ * @return the function value
+ * @throws FunctionEvaluationException if a runtime error occurs
+ * @see UnivariateRealFunction#value(double)
+ */
+ public double value(double z) throws FunctionEvaluationException {
+ try {
+ return evaluate(x, y, z);
+ } catch (DuplicateSampleAbscissaException e) {
+ throw new FunctionEvaluationException(z, e.getPattern(), e.getArguments(), e);
+ }
+ }
+
+ /**
+ * Returns the degree of the polynomial.
+ *
+ * @return the degree of the polynomial
+ */
+ public int degree() {
+ return x.length - 1;
+ }
+
+ /**
+ * Returns a copy of the interpolating points array.
+ * <p>
+ * Changes made to the returned copy will not affect the polynomial.
+ *
+ * @return a fresh copy of the interpolating points array
+ */
+ public double[] getInterpolatingPoints() {
+ double[] out = new double[x.length];
+ System.arraycopy(x, 0, out, 0, x.length);
+ return out;
+ }
+
+ /**
+ * Returns a copy of the interpolating values array.
+ * <p>
+ * Changes made to the returned copy will not affect the polynomial.
+ *
+ * @return a fresh copy of the interpolating values array
+ */
+ public double[] getInterpolatingValues() {
+ double[] out = new double[y.length];
+ System.arraycopy(y, 0, out, 0, y.length);
+ return out;
+ }
+
+ /**
+ * Returns a copy of the coefficients array.
+ * <p>
+ * Changes made to the returned copy will not affect the polynomial.
+ *
+ * @return a fresh copy of the coefficients array
+ */
+ public double[] getCoefficients() {
+ if (!coefficientsComputed) {
+ computeCoefficients();
+ }
+ double[] out = new double[coefficients.length];
+ System.arraycopy(coefficients, 0, out, 0, coefficients.length);
+ return out;
+ }
+
+ /**
+ * Evaluate the Lagrange polynomial using
+ * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
+ * Neville's Algorithm</a>. It takes O(N^2) time.
+ * <p>
+ * This function is made public static so that users can call it directly
+ * without instantiating PolynomialFunctionLagrangeForm object.
+ *
+ * @param x the interpolating points array
+ * @param y the interpolating values array
+ * @param z the point at which the function value is to be computed
+ * @return the function value
+ * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas
+ * @throws IllegalArgumentException if inputs are not valid
+ */
+ public static double evaluate(double x[], double y[], double z) throws
+ DuplicateSampleAbscissaException, IllegalArgumentException {
+
+ int i, j, n, nearest = 0;
+ double value, c[], d[], tc, td, divider, w, dist, min_dist;
+
+ verifyInterpolationArray(x, y);
+
+ n = x.length;
+ c = new double[n];
+ d = new double[n];
+ min_dist = Double.POSITIVE_INFINITY;
+ for (i = 0; i < n; i++) {
+ // initialize the difference arrays
+ c[i] = y[i];
+ d[i] = y[i];
+ // find out the abscissa closest to z
+ dist = Math.abs(z - x[i]);
+ if (dist < min_dist) {
+ nearest = i;
+ min_dist = dist;
+ }
+ }
+
+ // initial approximation to the function value at z
+ value = y[nearest];
+
+ for (i = 1; i < n; i++) {
+ for (j = 0; j < n-i; j++) {
+ tc = x[j] - z;
+ td = x[i+j] - z;
+ divider = x[j] - x[i+j];
+ if (divider == 0.0) {
+ // This happens only when two abscissas are identical.
+ throw new DuplicateSampleAbscissaException(x[i], i, i+j);
+ }
+ // update the difference arrays
+ w = (c[j+1] - d[j]) / divider;
+ c[j] = tc * w;
+ d[j] = td * w;
+ }
+ // sum up the difference terms to get the final value
+ if (nearest < 0.5*(n-i+1)) {
+ value += c[nearest]; // fork down
+ } else {
+ nearest--;
+ value += d[nearest]; // fork up
+ }
+ }
+
+ return value;
+ }
+
+ /**
+ * Calculate the coefficients of Lagrange polynomial from the
+ * interpolation data. It takes O(N^2) time.
+ * <p>
+ * Note this computation can be ill-conditioned. Use with caution
+ * and only when it is necessary.
+ *
+ * @throws ArithmeticException if any abscissas coincide
+ */
+ protected void computeCoefficients() throws ArithmeticException {
+ int i, j, n;
+ double c[], tc[], d, t;
+
+ n = degree() + 1;
+ coefficients = new double[n];
+ for (i = 0; i < n; i++) {
+ coefficients[i] = 0.0;
+ }
+
+ // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
+ c = new double[n+1];
+ c[0] = 1.0;
+ for (i = 0; i < n; i++) {
+ for (j = i; j > 0; j--) {
+ c[j] = c[j-1] - c[j] * x[i];
+ }
+ c[0] *= (-x[i]);
+ c[i+1] = 1;
+ }
+
+ tc = new double[n];
+ for (i = 0; i < n; i++) {
+ // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
+ d = 1;
+ for (j = 0; j < n; j++) {
+ if (i != j) {
+ d *= (x[i] - x[j]);
+ }
+ }
+ if (d == 0.0) {
+ // This happens only when two abscissas are identical.
+ throw new ArithmeticException
+ ("Identical abscissas cause division by zero.");
+ }
+ t = y[i] / d;
+ // Lagrange polynomial is the sum of n terms, each of which is a
+ // polynomial of degree n-1. tc[] are the coefficients of the i-th
+ // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
+ tc[n-1] = c[n]; // actually c[n] = 1
+ coefficients[n-1] += t * tc[n-1];
+ for (j = n-2; j >= 0; j--) {
+ tc[j] = c[j+1] + tc[j+1] * x[i];
+ coefficients[j] += t * tc[j];
+ }
+ }
+
+ coefficientsComputed = true;
+ }
+
+ /**
+ * Verifies that the interpolation arrays are valid.
+ * <p>
+ * The interpolating points must be distinct. However it is not
+ * verified here, it is checked in evaluate() and computeCoefficients().
+ *
+ * @throws IllegalArgumentException if not valid
+ * @see #evaluate(double[], double[], double)
+ * @see #computeCoefficients()
+ */
+ protected static void verifyInterpolationArray(double x[], double y[]) throws
+ IllegalArgumentException {
+
+ if (x.length < 2 || y.length < 2) {
+ throw new IllegalArgumentException
+ ("Interpolation requires at least two points.");
+ }
+ if (x.length != y.length) {
+ throw new IllegalArgumentException
+ ("Abscissa and value arrays must have the same length.");
+ }
+ }
+}
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
('svn:executable' removed)
Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java?view=auto&rev=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java Mon Feb 12 11:27:16 2007
@@ -0,0 +1,38 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math;
+
+import java.util.Locale;
+
+import junit.framework.TestCase;
+
+/**
+ * @version $Revision:$
+ */
+public class DuplicateSampleAbscissaExceptionTest extends TestCase {
+
+ public void testConstructor(){
+ DuplicateSampleAbscissaException ex = new DuplicateSampleAbscissaException(1.2, 10, 11);
+ assertNull(ex.getCause());
+ assertNotNull(ex.getMessage());
+ assertTrue(ex.getMessage().indexOf("1.2") > 0);
+ assertEquals(1.2, ex.getDuplicateAbscissa(), 0);
+ assertFalse(ex.getMessage().equals(ex.getMessage(Locale.FRENCH)));
+ }
+
+}
Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
------------------------------------------------------------------------------
svn:eol-style = native
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org