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 }