View Javadoc

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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 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          double a[], c[];
61  
62          PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
63  
64          /**
65           * When used for interpolation, the Newton form formula becomes
66           * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
67           *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
68           * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
69           * <p>
70           * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
71           */
72          c = new double[x.length-1];
73          for (int i = 0; i < c.length; i++) {
74              c[i] = x[i];
75          }
76          a = computeDividedDifference(x, y);
77  
78          return new PolynomialFunctionNewtonForm(a, c);
79  
80      }
81  
82      /**
83       * Returns a copy of the divided difference array.
84       * <p> 
85       * The divided difference array is defined recursively by <pre>
86       * f[x0] = f(x0)
87       * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
88       * </pre></p>
89       * <p>
90       * The computational complexity is O(N^2).</p>
91       *
92       * @param x the interpolating points array
93       * @param y the interpolating values array
94       * @return a fresh copy of the divided difference array
95       * @throws DuplicateSampleAbscissaException if any abscissas coincide
96       */
97      protected static double[] computeDividedDifference(double x[], double y[])
98          throws DuplicateSampleAbscissaException {
99  
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 }