Home | History | Annotate | Download | only in interpolation 
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) under one or more
      3  * contributor license agreements.  See the NOTICE file distributed with
      4  * this work for additional information regarding copyright ownership.
      5  * The ASF licenses this file to You under the Apache License, Version 2.0
      6  * (the "License"); you may not use this file except in compliance with
      7  * the License.  You may obtain a copy of the License at
      8  *
      9  *      http://www.apache.org/licenses/LICENSE-2.0
     10  *
     11  * Unless required by applicable law or agreed to in writing, software
     12  * distributed under the License is distributed on an "AS IS" BASIS,
     13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     14  * See the License for the specific language governing permissions and
     15  * limitations under the License.
     16  */
     17 package org.apache.commons.math.analysis.interpolation;
     18 
     19 import java.io.Serializable;
     20 
     21 import org.apache.commons.math.DuplicateSampleAbscissaException;
     22 import org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm;
     23 import org.apache.commons.math.analysis.polynomials.PolynomialFunctionNewtonForm;
     24 
     25 /**
     26  * Implements the <a href="
     27  * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
     28  * Divided Difference Algorithm</a> for interpolation of real univariate
     29  * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
     30  * ISBN 038795452X, chapter 2.
     31  * <p>
     32  * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
     33  * this class provides an easy-to-use interface to it.</p>
     34  *
     35  * @version $Revision: 825919 $ $Date: 2009-10-16 16:51:55 +0200 (ven. 16 oct. 2009) $
     36  * @since 1.2
     37  */
     38 public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
     39     Serializable {
     40 
     41     /** serializable version identifier */
     42     private static final long serialVersionUID = 107049519551235069L;
     43 
     44     /**
     45      * Computes an interpolating function for the data set.
     46      *
     47      * @param x the interpolating points array
     48      * @param y the interpolating values array
     49      * @return a function which interpolates the data set
     50      * @throws DuplicateSampleAbscissaException if arguments are invalid
     51      */
     52     public PolynomialFunctionNewtonForm interpolate(double x[], double y[]) throws
     53         DuplicateSampleAbscissaException {
     54 
     55         /**
     56          * a[] and c[] are defined in the general formula of Newton form:
     57          * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
     58          *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
     59          */
     60         PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
     61 
     62         /**
     63          * When used for interpolation, the Newton form formula becomes
     64          * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
     65          *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
     66          * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
     67          * <p>
     68          * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
     69          */
     70         final double[] c = new double[x.length-1];
     71         System.arraycopy(x, 0, c, 0, c.length);
     72 
     73         final double[] a = computeDividedDifference(x, y);
     74         return new PolynomialFunctionNewtonForm(a, c);
     75 
     76     }
     77 
     78     /**
     79      * Returns a copy of the divided difference array.
     80      * <p>
     81      * The divided difference array is defined recursively by <pre>
     82      * f[x0] = f(x0)
     83      * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
     84      * </pre></p>
     85      * <p>
     86      * The computational complexity is O(N^2).</p>
     87      *
     88      * @param x the interpolating points array
     89      * @param y the interpolating values array
     90      * @return a fresh copy of the divided difference array
     91      * @throws DuplicateSampleAbscissaException if any abscissas coincide
     92      */
     93     protected static double[] computeDividedDifference(final double x[], final double y[])
     94         throws DuplicateSampleAbscissaException {
     95 
     96         PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
     97 
     98         final double[] divdiff = y.clone(); // initialization
     99 
    100         final int n = x.length;
    101         final double[] a = new double [n];
    102         a[0] = divdiff[0];
    103         for (int i = 1; i < n; i++) {
    104             for (int j = 0; j < n-i; j++) {
    105                 final double denominator = x[j+i] - x[j];
    106                 if (denominator == 0.0) {
    107                     // This happens only when two abscissas are identical.
    108                     throw new DuplicateSampleAbscissaException(x[j], j, j+i);
    109                 }
    110                 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
    111             }
    112             a[i] = divdiff[0];
    113         }
    114 
    115         return a;
    116     }
    117 }
    118