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.math.analysis.interpolation;
018    
019    import java.io.Serializable;
020    
021    import org.apache.commons.math.DuplicateSampleAbscissaException;
022    import org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm;
023    import org.apache.commons.math.analysis.polynomials.PolynomialFunctionNewtonForm;
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     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
036     * @since 1.2
037     */
038    public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
039        Serializable {
040    
041        /** serializable version identifier */
042        private static final long serialVersionUID = 107049519551235069L;
043    
044        /**
045         * Computes an interpolating function for the data set.
046         *
047         * @param x the interpolating points array
048         * @param y the interpolating values array
049         * @return a function which interpolates the data set
050         * @throws DuplicateSampleAbscissaException if arguments are invalid
051         */
052        public PolynomialFunctionNewtonForm interpolate(double x[], double y[]) throws
053            DuplicateSampleAbscissaException {
054    
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            double a[], c[];
061    
062            PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
063    
064            /**
065             * When used for interpolation, the Newton form formula becomes
066             * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
067             *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
068             * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
069             * <p>
070             * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
071             */
072            c = new double[x.length-1];
073            for (int i = 0; i < c.length; i++) {
074                c[i] = x[i];
075            }
076            a = computeDividedDifference(x, y);
077    
078            return new PolynomialFunctionNewtonForm(a, c);
079    
080        }
081    
082        /**
083         * Returns a copy of the divided difference array.
084         * <p> 
085         * The divided difference array is defined recursively by <pre>
086         * f[x0] = f(x0)
087         * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
088         * </pre></p>
089         * <p>
090         * The computational complexity is O(N^2).</p>
091         *
092         * @param x the interpolating points array
093         * @param y the interpolating values array
094         * @return a fresh copy of the divided difference array
095         * @throws DuplicateSampleAbscissaException if any abscissas coincide
096         */
097        protected static double[] computeDividedDifference(double x[], double y[])
098            throws DuplicateSampleAbscissaException {
099    
100            int i, j, n;
101            double divdiff[], a[], denominator;
102    
103            PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
104    
105            n = x.length;
106            divdiff = new double[n];
107            for (i = 0; i < n; i++) {
108                divdiff[i] = y[i];      // initialization
109            }
110    
111            a = new double [n];
112            a[0] = divdiff[0];
113            for (i = 1; i < n; i++) {
114                for (j = 0; j < n-i; j++) {
115                    denominator = x[j+i] - x[j];
116                    if (denominator == 0.0) {
117                        // This happens only when two abscissas are identical.
118                        throw new DuplicateSampleAbscissaException(x[j], j, j+i);
119                    }
120                    divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
121                }
122                a[i] = divdiff[0];
123            }
124    
125            return a;
126        }
127    }