/*
* 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.polynomials;
import org.apache.commons.math.DuplicateSampleAbscissaException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* 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.</p>
*
* @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 févr. 2011) $
* @since 1.2
*/
public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
/**
* 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).
*/
private final double x[];
/**
* Function values at interpolating points.
*/
private final double 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.</p>
*
* @param x interpolating points
* @param y function values at interpolating points
* @throws IllegalArgumentException if input arrays are not valid
*/
public 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;
}
/** {@inheritDoc} */
public double value(double z) throws FunctionEvaluationException {
try {
return evaluate(x, y, z);
} catch (DuplicateSampleAbscissaException e) {
throw new FunctionEvaluationException(z, e.getSpecificPattern(), e.getGeneralPattern(), e.getArguments());
}
}
/**
* 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.</p>
*
* @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.</p>
*
* @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.</p>
* <p>
* Note that coefficients computation can be ill-conditioned. Use with caution
* and only when it is necessary.</p>
*
* @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.</p>
*
* @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 {
verifyInterpolationArray(x, y);
int nearest = 0;
final int n = x.length;
final double[] c = new double[n];
final double[] d = new double[n];
double min_dist = Double.POSITIVE_INFINITY;
for (int i = 0; i < n; i++) {
// initialize the difference arrays
c[i] = y[i];
d[i] = y[i];
// find out the abscissa closest to z
final double dist = FastMath.abs(z - x[i]);
if (dist < min_dist) {
nearest = i;
min_dist = dist;
}
}
// initial approximation to the function value at z
double value = y[nearest];
for (int i = 1; i < n; i++) {
for (int j = 0; j < n-i; j++) {
final double tc = x[j] - z;
final double td = x[i+j] - z;
final double 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
final double 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.</p>
*
* @throws ArithmeticException if any abscissas coincide
*/
protected void computeCoefficients() throws ArithmeticException {
final int n = degree() + 1;
coefficients = new double[n];
for (int 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])
final double[] c = new double[n+1];
c[0] = 1.0;
for (int i = 0; i < n; i++) {
for (int j = i; j > 0; j--) {
c[j] = c[j-1] - c[j] * x[i];
}
c[0] *= -x[i];
c[i+1] = 1;
}
final double[] tc = new double[n];
for (int 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])
double d = 1;
for (int 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.
for (int k = 0; k < n; ++k) {
if ((i != k) && (x[i] == x[k])) {
throw MathRuntimeException.createArithmeticException(
LocalizedFormats.IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO,
i, k, x[i]);
}
}
}
final double 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 (int 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 arrays features checked by this method are that both arrays have the
* same length and this length is at least 2.
* </p>
* <p>
* The interpolating points must be distinct. However it is not
* verified here, it is checked in evaluate() and computeCoefficients().
* </p>
*
* @param x the interpolating points array
* @param y the interpolating values array
* @throws IllegalArgumentException if not valid
* @see #evaluate(double[], double[], double)
* @see #computeCoefficients()
*/
public static void verifyInterpolationArray(double x[], double y[])
throws IllegalArgumentException {
if (x.length != y.length) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, x.length, y.length);
}
if (x.length < 2) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length);
}
}
}