    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.interpolation;
    018
    019    import java.io.Serializable;
    020    import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm;
    021    import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionNewtonForm;
    022    import org.apache.commons.math3.exception.DimensionMismatchException;
    023    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    024    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
    025
    026    /**
    027     * Implements the <a href="
    028     * http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
    029     * Divided Difference Algorithm</a> for interpolation of real univariate
    030     * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
    031     * ISBN 038795452X, chapter 2.
    032     * <p>
    033     * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
    034     * this class provides an easy-to-use interface to it.</p>
    035     *
    036     * @version $Id: DividedDifferenceInterpolator.java 1385313 2012-09-16 16:35:23Z tn $
    037     * @since 1.2
    038     */
    039    public class DividedDifferenceInterpolator
    040        implements UnivariateInterpolator, Serializable {
    041        /** serializable version identifier */
    042        private static final long serialVersionUID = 107049519551235069L;
    043
    044        /**
    045         * Compute an interpolating function for the dataset.
    046         *
    047         * @param x Interpolating points array.
    048         * @param y Interpolating values array.
    049         * @return a function which interpolates the dataset.
    050         * @throws DimensionMismatchException if the array lengths are different.
    051         * @throws NumberIsTooSmallException if the number of points is less than 2.
    052         * @throws NonMonotonicSequenceException if {@code x} is not sorted in
    053         * strictly increasing order.
    054         */
    055        public PolynomialFunctionNewtonForm interpolate(double x[], double y[])
    056            throws DimensionMismatchException,
    057                   NumberIsTooSmallException,
    058                   NonMonotonicSequenceException {
    059            /**
    060             * a[] and c[] are defined in the general formula of Newton form:
    061             * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
    062             *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
    063             */
    064            PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
    065
    066            /**
    067             * When used for interpolation, the Newton form formula becomes
    068             * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
    069             *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
    070             * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
    071             * <p>
    072             * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
    073             */
    074            final double[] c = new double[x.length-1];
    075            System.arraycopy(x, 0, c, 0, c.length);
    076
    077            final double[] a = computeDividedDifference(x, y);
    078            return new PolynomialFunctionNewtonForm(a, c);
    079        }
    080
    081        /**
    082         * Return a copy of the divided difference array.
    083         * <p>
    084         * The divided difference array is defined recursively by <pre>
    085         * f[x0] = f(x0)
    086         * f[x0,x1,...,xk] = (f[x1,...,xk] - f[x0,...,x[k-1]]) / (xk - x0)
    087         * </pre></p>
    088         * <p>
    089         * The computational complexity is O(N^2).</p>
    090         *
    091         * @param x Interpolating points array.
    092         * @param y Interpolating values array.
    093         * @return a fresh copy of the divided difference array.
    094         * @throws DimensionMismatchException if the array lengths are different.
    095         * @throws NumberIsTooSmallException if the number of points is less than 2.
    096         * @throws NonMonotonicSequenceException
    097         * if {@code x} is not sorted in strictly increasing order.
    098         */
    099        protected static double[] computeDividedDifference(final double x[], final double y[])
    100            throws DimensionMismatchException,
    101                   NumberIsTooSmallException,
    102                   NonMonotonicSequenceException {
    103            PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
    104
    105            final double[] divdiff = y.clone(); // initialization
    106
    107            final int n = x.length;
    108            final double[] a = new double [n];
    109            a[0] = divdiff[0];
    110            for (int i = 1; i < n; i++) {
    111                for (int j = 0; j < n-i; j++) {
    112                    final double denominator = x[j+i] - x[j];
    113                    divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
    114                }
    115                a[i] = divdiff[0];
    116            }
    117
    118            return a;
    119        }
    120    }