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 import org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm;
021 import org.apache.commons.math.analysis.polynomials.PolynomialFunctionNewtonForm;
022
023 /**
024 * Implements the <a href="
025 * http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
026 * Divided Difference Algorithm</a> for interpolation of real univariate
027 * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
028 * ISBN 038795452X, chapter 2.
029 * <p>
030 * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
031 * this class provides an easy-to-use interface to it.</p>
032 *
033 * @version $Id: DividedDifferenceInterpolator.java 1179928 2011-10-07 03:20:39Z psteitz $
034 * @since 1.2
035 */
036 public class DividedDifferenceInterpolator
037 implements UnivariateRealInterpolator, Serializable {
038 /** serializable version identifier */
039 private static final long serialVersionUID = 107049519551235069L;
040
041 /**
042 * Compute an interpolating function for the dataset.
043 *
044 * @param x Interpolating points array.
045 * @param y Interpolating values array.
046 * @return a function which interpolates the dataset.
047 * @throws org.apache.commons.math.exception.DimensionMismatchException
048 * if the array lengths are different.
049 * @throws org.apache.commons.math.exception.NumberIsTooSmallException
050 * if the number of points is less than 2.
051 * @throws org.apache.commons.math.exception.NonMonotonicSequenceException
052 * if {@code x} is not sorted in strictly increasing order.
053 */
054 public PolynomialFunctionNewtonForm interpolate(double x[], double y[]) {
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></p>
084 * <p>
085 * The computational complexity is O(N^2).</p>
086 *
087 * @param x Interpolating points array.
088 * @param y Interpolating values array.
089 * @return a fresh copy of the divided difference array.
090 * @throws org.apache.commons.math.exception.DimensionMismatchException
091 * if the array lengths are different.
092 * @throws org.apache.commons.math.exception.NumberIsTooSmallException
093 * if the number of points is less than 2.
094 * @throws org.apache.commons.math.exception.NonMonotonicSequenceException
095 * if {@code x} is not sorted in strictly increasing order.
096 */
097 protected static double[] computeDividedDifference(final double x[], final double y[]) {
098 PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
099
100 final double[] divdiff = y.clone(); // initialization
101
102 final int n = x.length;
103 final double[] a = new double [n];
104 a[0] = divdiff[0];
105 for (int i = 1; i < n; i++) {
106 for (int j = 0; j < n-i; j++) {
107 final double denominator = x[j+i] - x[j];
108 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
109 }
110 a[i] = divdiff[0];
111 }
112
113 return a;
114 }
115 }