Coverage Report - org.apache.commons.math.analysis.DividedDifferenceInterpolator

Classes in this File Line Coverage Branch Coverage Complexity
DividedDifferenceInterpolator
100% 
100% 
4

 1  
 /*
 2  
  * Copyright 2003-2005 The Apache Software Foundation.
 3  
  *
 4  
  * Licensed under the Apache License, Version 2.0 (the "License");
 5  
  * you may not use this file except in compliance with the License.
 6  
  * You may obtain a copy of the License at
 7  
  *
 8  
  *      http://www.apache.org/licenses/LICENSE-2.0
 9  
  *
 10  
  * Unless required by applicable law or agreed to in writing, software
 11  
  * distributed under the License is distributed on an "AS IS" BASIS,
 12  
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 13  
  * See the License for the specific language governing permissions and
 14  
  * limitations under the License.
 15  
  */
 16  
 package org.apache.commons.math.analysis;
 17  
 
 18  
 import java.io.Serializable;
 19  
 import org.apache.commons.math.MathException;
 20  
 
 21  
 /**
 22  
  * Implements the <a href="
 23  
  * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
 24  
  * Divided Difference Algorithm</a> for interpolation of real univariate
 25  
  * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
 26  
  * ISBN 038795452X, chapter 2.
 27  
  * <p>
 28  
  * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
 29  
  * this class provides an easy-to-use interface to it.
 30  
  *
 31  
  * @version $Revision$ $Date$
 32  
  */
 33  6
 public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
 34  
     Serializable {
 35  
 
 36  
     /** serializable version identifier */
 37  
     static final long serialVersionUID = 107049519551235069L;
 38  
 
 39  
     /**
 40  
      * Computes an interpolating function for the data set.
 41  
      *
 42  
      * @param x the interpolating points array
 43  
      * @param y the interpolating values array
 44  
      * @return a function which interpolates the data set
 45  
      * @throws MathException if arguments are invalid
 46  
      */
 47  
     public UnivariateRealFunction interpolate(double x[], double y[]) throws
 48  
         MathException {
 49  
 
 50  
         /**
 51  
          * a[] and c[] are defined in the general formula of Newton form:
 52  
          * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
 53  
          *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
 54  
          */
 55  
         double a[], c[];
 56  
 
 57  6
         PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
 58  
 
 59  
         /**
 60  
          * When used for interpolation, the Newton form formula becomes
 61  
          * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
 62  
          *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
 63  
          * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
 64  
          * <p>
 65  
          * Note x[], y[], a[] have the same length but c[]'s size is one less.
 66  
          */
 67  6
         c = new double[x.length-1];
 68  30
         for (int i = 0; i < c.length; i++) {
 69  24
             c[i] = x[i];
 70  
         }
 71  6
         a = computeDividedDifference(x, y);
 72  
 
 73  
         PolynomialFunctionNewtonForm p;
 74  4
         p = new PolynomialFunctionNewtonForm(a, c);
 75  4
         return p;
 76  
     }
 77  
 
 78  
     /**
 79  
      * Returns a copy of the divided difference array.
 80  
      * <p> 
 81  
      * The divided difference array is defined recursively by <pre>
 82  
      * f[x0] = f(x0)
 83  
      * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
 84  
      * </pre><p>
 85  
      * The computational complexity is O(N^2).
 86  
      *
 87  
      * @return a fresh copy of the divided difference array
 88  
      * @throws MathException if any abscissas coincide
 89  
      */
 90  
     protected static double[] computeDividedDifference(double x[], double y[])
 91  
         throws MathException {
 92  
 
 93  
         int i, j, n;
 94  
         double divdiff[], a[], denominator;
 95  
 
 96  6
         PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
 97  
 
 98  6
         n = x.length;
 99  6
         divdiff = new double[n];
 100  36
         for (i = 0; i < n; i++) {
 101  30
             divdiff[i] = y[i];      // initialization
 102  
         }
 103  
 
 104  6
         a = new double [n];
 105  6
         a[0] = divdiff[0];
 106  24
         for (i = 1; i < n; i++) {
 107  72
             for (j = 0; j < n-i; j++) {
 108  54
                 denominator = x[j+i] - x[j];
 109  54
                 if (denominator == 0.0) {
 110  
                     // This happens only when two abscissas are identical.
 111  2
                     throw new MathException
 112  
                         ("Identical abscissas cause division by zero.");
 113  
                 }
 114  52
                 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
 115  
             }
 116  18
             a[i] = divdiff[0];
 117  
         }
 118  
 
 119  4
         return a;
 120  
     }
 121  
 }