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; 18 19 import java.io.Serializable; 20 21 import org.apache.commons.math.DuplicateSampleAbscissaException; 22 23 /** 24 * Implements the <a href=" 25 * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html"> 26 * Divided Difference Algorithm</a> for interpolation of real univariate 27 * functions. For reference, see <b>Introduction to Numerical Analysis</b>, 28 * ISBN 038795452X, chapter 2. 29 * <p> 30 * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm, 31 * this class provides an easy-to-use interface to it.</p> 32 * 33 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 34 * @since 1.2 35 */ 36 public class DividedDifferenceInterpolator implements UnivariateRealInterpolator, 37 Serializable { 38 39 /** serializable version identifier */ 40 private static final long serialVersionUID = 107049519551235069L; 41 42 /** 43 * Computes an interpolating function for the data set. 44 * 45 * @param x the interpolating points array 46 * @param y the interpolating values array 47 * @return a function which interpolates the data set 48 * @throws DuplicateSampleAbscissaException if arguments are invalid 49 */ 50 public UnivariateRealFunction interpolate(double x[], double y[]) throws 51 DuplicateSampleAbscissaException { 52 53 /** 54 * a[] and c[] are defined in the general formula of Newton form: 55 * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + 56 * a[n](x-c[0])(x-c[1])...(x-c[n-1]) 57 */ 58 double a[], c[]; 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 c = new double[x.length-1]; 71 for (int i = 0; i < c.length; i++) { 72 c[i] = x[i]; 73 } 74 a = computeDividedDifference(x, y); 75 76 PolynomialFunctionNewtonForm p; 77 p = new PolynomialFunctionNewtonForm(a, c); 78 return p; 79 } 80 81 /** 82 * Returns a copy of the divided difference array. 83 * <p> 84 * The divided difference array is defined recursively by <pre> 85 * f[x0] = f(x0) 86 * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0) 87 * </pre></p> 88 * <p> 89 * The computational complexity is O(N^2).</p> 90 * 91 * @param x the interpolating points array 92 * @param y the interpolating values array 93 * @return a fresh copy of the divided difference array 94 * @throws DuplicateSampleAbscissaException if any abscissas coincide 95 */ 96 protected static double[] computeDividedDifference(double x[], double y[]) 97 throws DuplicateSampleAbscissaException { 98 99 int i, j, n; 100 double divdiff[], a[], denominator; 101 102 PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y); 103 104 n = x.length; 105 divdiff = new double[n]; 106 for (i = 0; i < n; i++) { 107 divdiff[i] = y[i]; // initialization 108 } 109 110 a = new double [n]; 111 a[0] = divdiff[0]; 112 for (i = 1; i < n; i++) { 113 for (j = 0; j < n-i; j++) { 114 denominator = x[j+i] - x[j]; 115 if (denominator == 0.0) { 116 // This happens only when two abscissas are identical. 117 throw new DuplicateSampleAbscissaException(x[j], j, j+i); 118 } 119 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator; 120 } 121 a[i] = divdiff[0]; 122 } 123 124 return a; 125 } 126 }