Added and used a specialized exception for duplicate abscissas in sampled functions
git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@506592 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
21a95478c2
commit
1bd2978235
|
|
@ -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();
|
||||
}
|
||||
|
||||
}
|
||||
244
src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
Executable file → Normal file
244
src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
Executable file → Normal file
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
586
src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
Executable file → Normal file
586
src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
Executable file → Normal file
|
|
@ -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.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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)));
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue