diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java')
-rw-r--r-- | src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java | 239 |
1 files changed, 239 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java new file mode 100644 index 0000000..15ed322 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java @@ -0,0 +1,239 @@ +/* + * 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.math3.analysis.interpolation; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.exception.MathArithmeticException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.CombinatoricsUtils; + +/** Polynomial interpolator using both sample values and sample derivatives. + * <p> + * The interpolation polynomials match all sample points, including both values + * and provided derivatives. There is one polynomial for each component of + * the values vector. All polynomials have the same degree. The degree of the + * polynomials depends on the number of points and number of derivatives at each + * point. For example the interpolation polynomials for n sample points without + * any derivatives all have degree n-1. The interpolation polynomials for n + * sample points with the two extreme points having value and first derivative + * and the remaining points having value only all have degree n+1. The + * interpolation polynomial for n sample points with value, first and second + * derivative for all points all have degree 3n-1. + * </p> + * + * @since 3.1 + */ +public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction { + + /** Sample abscissae. */ + private final List<Double> abscissae; + + /** Top diagonal of the divided differences array. */ + private final List<double[]> topDiagonal; + + /** Bottom diagonal of the divided differences array. */ + private final List<double[]> bottomDiagonal; + + /** Create an empty interpolator. + */ + public HermiteInterpolator() { + this.abscissae = new ArrayList<Double>(); + this.topDiagonal = new ArrayList<double[]>(); + this.bottomDiagonal = new ArrayList<double[]>(); + } + + /** Add a sample point. + * <p> + * This method must be called once for each sample point. It is allowed to + * mix some calls with values only with calls with values and first + * derivatives. + * </p> + * <p> + * The point abscissae for all calls <em>must</em> be different. + * </p> + * @param x abscissa of the sample point + * @param value value and derivatives of the sample point + * (if only one row is passed, it is the value, if two rows are + * passed the first one is the value and the second the derivative + * and so on) + * @exception ZeroException if the abscissa difference between added point + * and a previous point is zero (i.e. the two points are at same abscissa) + * @exception MathArithmeticException if the number of derivatives is larger + * than 20, which prevents computation of a factorial + */ + public void addSamplePoint(final double x, final double[] ... value) + throws ZeroException, MathArithmeticException { + + for (int i = 0; i < value.length; ++i) { + + final double[] y = value[i].clone(); + if (i > 1) { + double inv = 1.0 / CombinatoricsUtils.factorial(i); + for (int j = 0; j < y.length; ++j) { + y[j] *= inv; + } + } + + // update the bottom diagonal of the divided differences array + final int n = abscissae.size(); + bottomDiagonal.add(n - i, y); + double[] bottom0 = y; + for (int j = i; j < n; ++j) { + final double[] bottom1 = bottomDiagonal.get(n - (j + 1)); + final double inv = 1.0 / (x - abscissae.get(n - (j + 1))); + if (Double.isInfinite(inv)) { + throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x); + } + for (int k = 0; k < y.length; ++k) { + bottom1[k] = inv * (bottom0[k] - bottom1[k]); + } + bottom0 = bottom1; + } + + // update the top diagonal of the divided differences array + topDiagonal.add(bottom0.clone()); + + // update the abscissae array + abscissae.add(x); + + } + + } + + /** Compute the interpolation polynomials. + * @return interpolation polynomials array + * @exception NoDataException if sample is empty + */ + public PolynomialFunction[] getPolynomials() + throws NoDataException { + + // safety check + checkInterpolation(); + + // iteration initialization + final PolynomialFunction zero = polynomial(0); + PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length]; + for (int i = 0; i < polynomials.length; ++i) { + polynomials[i] = zero; + } + PolynomialFunction coeff = polynomial(1); + + // build the polynomials by iterating on the top diagonal of the divided differences array + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] tdi = topDiagonal.get(i); + for (int k = 0; k < polynomials.length; ++k) { + polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k]))); + } + coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0)); + } + + return polynomials; + + } + + /** Interpolate value at a specified abscissa. + * <p> + * Calling this method is equivalent to call the {@link PolynomialFunction#value(double) + * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials}, + * except it does not build the intermediate polynomials, so this method is faster and + * numerically more stable. + * </p> + * @param x interpolation abscissa + * @return interpolated value + * @exception NoDataException if sample is empty + */ + public double[] value(double x) + throws NoDataException { + + // safety check + checkInterpolation(); + + final double[] value = new double[topDiagonal.get(0).length]; + double valueCoeff = 1; + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < value.length; ++k) { + value[k] += dividedDifference[k] * valueCoeff; + } + final double deltaX = x - abscissae.get(i); + valueCoeff *= deltaX; + } + + return value; + + } + + /** Interpolate value at a specified abscissa. + * <p> + * Calling this method is equivalent to call the {@link + * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials + * returned by {@link #getPolynomials() getPolynomials}, except it does not build the + * intermediate polynomials, so this method is faster and numerically more stable. + * </p> + * @param x interpolation abscissa + * @return interpolated value + * @exception NoDataException if sample is empty + */ + public DerivativeStructure[] value(final DerivativeStructure x) + throws NoDataException { + + // safety check + checkInterpolation(); + + final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length]; + Arrays.fill(value, x.getField().getZero()); + DerivativeStructure valueCoeff = x.getField().getOne(); + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < value.length; ++k) { + value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k])); + } + final DerivativeStructure deltaX = x.subtract(abscissae.get(i)); + valueCoeff = valueCoeff.multiply(deltaX); + } + + return value; + + } + + /** Check interpolation can be performed. + * @exception NoDataException if interpolation cannot be performed + * because sample is empty + */ + private void checkInterpolation() throws NoDataException { + if (abscissae.isEmpty()) { + throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE); + } + } + + /** Create a polynomial from its coefficients. + * @param c polynomials coefficients + * @return polynomial + */ + private PolynomialFunction polynomial(double ... c) { + return new PolynomialFunction(c); + } + +} |