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 }