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 */
017package org.apache.commons.math4.legacy.analysis.interpolation;
018
019import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunctionLagrangeForm;
020import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunctionNewtonForm;
021import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
022import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
023import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
024
025/**
026 * Implements the <a href=
027 * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
028 * Divided Difference Algorithm</a> for interpolation of real univariate
029 * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
030 * ISBN 038795452X, chapter 2.
031 * <p>
032 * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
033 * this class provides an easy-to-use interface to it.</p>
034 *
035 * @since 1.2
036 */
037public class DividedDifferenceInterpolator
038    implements UnivariateInterpolator {
039    /**
040     * Compute an interpolating function for the dataset.
041     *
042     * @param x Interpolating points array.
043     * @param y Interpolating values array.
044     * @return a function which interpolates the dataset.
045     * @throws DimensionMismatchException if the array lengths are different.
046     * @throws NumberIsTooSmallException if the number of points is less than 2.
047     * @throws NonMonotonicSequenceException if {@code x} is not sorted in
048     * strictly increasing order.
049     */
050    @Override
051    public PolynomialFunctionNewtonForm interpolate(double[] x, double[] y)
052        throws DimensionMismatchException,
053               NumberIsTooSmallException,
054               NonMonotonicSequenceException {
055        /*
056         * a[] and c[] are defined in the general formula of Newton form:
057         * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
058         *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
059         */
060        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
061
062        /*
063         * When used for interpolation, the Newton form formula becomes
064         * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
065         *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
066         * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
067         * <p>
068         * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
069         */
070        final double[] c = new double[x.length-1];
071        System.arraycopy(x, 0, c, 0, c.length);
072
073        final double[] a = computeDividedDifference(x, y);
074        return new PolynomialFunctionNewtonForm(a, c);
075    }
076
077    /**
078     * Return a copy of the divided difference array.
079     * <p>
080     * The divided difference array is defined recursively by <pre>
081     * f[x0] = f(x0)
082     * f[x0,x1,...,xk] = (f[x1,...,xk] - f[x0,...,x[k-1]]) / (xk - x0)
083     * </pre>
084     * <p>
085     * The computational complexity is \(O(n^2)\) where \(n\) is the common
086     * length of {@code x} and {@code y}.</p>
087     *
088     * @param x Interpolating points array.
089     * @param y Interpolating values array.
090     * @return a fresh copy of the divided difference array.
091     * @throws DimensionMismatchException if the array lengths are different.
092     * @throws NumberIsTooSmallException if the number of points is less than 2.
093     * @throws NonMonotonicSequenceException
094     * if {@code x} is not sorted in strictly increasing order.
095     */
096    protected static double[] computeDividedDifference(final double[] x, final double[] y)
097        throws DimensionMismatchException,
098               NumberIsTooSmallException,
099               NonMonotonicSequenceException {
100        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
101
102        final double[] divdiff = y.clone(); // initialization
103
104        final int n = x.length;
105        final double[] a = new double [n];
106        a[0] = divdiff[0];
107        for (int i = 1; i < n; i++) {
108            for (int j = 0; j < n-i; j++) {
109                final double denominator = x[j+i] - x[j];
110                divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
111            }
112            a[i] = divdiff[0];
113        }
114
115        return a;
116    }
117}