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

Classes in this File Line Coverage Branch Coverage Complexity
PolynomialFunctionLagrangeForm
91% 
100% 
3.333

 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.FunctionEvaluationException;
 20  
 
 21  
 /**
 22  
  * Implements the representation of a real polynomial function in
 23  
  * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
 24  
  * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
 25  
  * Analysis</b>, ISBN 038795452X, chapter 2.
 26  
  * <p>
 27  
  * The approximated function should be smooth enough for Lagrange polynomial
 28  
  * to work well. Otherwise, consider using splines instead.
 29  
  *
 30  
  * @version $Revision$ $Date$
 31  
  */
 32  
 public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
 33  
     Serializable {
 34  
 
 35  
     /** serializable version identifier */
 36  
     static final long serialVersionUID = -3965199246151093920L;
 37  
 
 38  
     /**
 39  
      * The coefficients of the polynomial, ordered by degree -- i.e.
 40  
      * coefficients[0] is the constant term and coefficients[n] is the 
 41  
      * coefficient of x^n where n is the degree of the polynomial.
 42  
      */
 43  
     private double coefficients[];
 44  
 
 45  
     /**
 46  
      * Interpolating points (abscissas) and the function values at these points.
 47  
      */
 48  
     private double x[], y[];
 49  
 
 50  
     /**
 51  
      * Whether the polynomial coefficients are available.
 52  
      */
 53  
     private boolean coefficientsComputed;
 54  
 
 55  
     /**
 56  
      * Construct a Lagrange polynomial with the given abscissas and function
 57  
      * values. The order of interpolating points are not important.
 58  
      * <p>
 59  
      * The constructor makes copy of the input arrays and assigns them.
 60  
      * 
 61  
      * @param x interpolating points
 62  
      * @param y function values at interpolating points
 63  
      * @throws IllegalArgumentException if input arrays are not valid
 64  
      */
 65  
     PolynomialFunctionLagrangeForm(double x[], double y[]) throws
 66  16
         IllegalArgumentException {
 67  
 
 68  16
         verifyInterpolationArray(x, y);
 69  12
         this.x = new double[x.length];
 70  12
         this.y = new double[y.length];
 71  12
         System.arraycopy(x, 0, this.x, 0, x.length);
 72  12
         System.arraycopy(y, 0, this.y, 0, y.length);
 73  12
         coefficientsComputed = false;
 74  12
     }
 75  
 
 76  
     /**
 77  
      * Calculate the function value at the given point.
 78  
      *
 79  
      * @param z the point at which the function value is to be computed
 80  
      * @return the function value
 81  
      * @throws FunctionEvaluationException if a runtime error occurs
 82  
      * @see UnivariateRealFunction#value(double)
 83  
      */
 84  
     public double value(double z) throws FunctionEvaluationException {
 85  30
        return evaluate(x, y, z);
 86  
     }
 87  
 
 88  
     /**
 89  
      * Returns the degree of the polynomial.
 90  
      * 
 91  
      * @return the degree of the polynomial
 92  
      */
 93  
     public int degree() {
 94  12
         return x.length - 1;
 95  
     }
 96  
 
 97  
     /**
 98  
      * Returns a copy of the interpolating points array.
 99  
      * <p>
 100  
      * Changes made to the returned copy will not affect the polynomial.
 101  
      * 
 102  
      * @return a fresh copy of the interpolating points array
 103  
      */
 104  
     public double[] getInterpolatingPoints() {
 105  0
         double[] out = new double[x.length];
 106  0
         System.arraycopy(x, 0, out, 0, x.length);
 107  0
         return out;
 108  
     }
 109  
 
 110  
     /**
 111  
      * Returns a copy of the interpolating values array.
 112  
      * <p>
 113  
      * Changes made to the returned copy will not affect the polynomial.
 114  
      * 
 115  
      * @return a fresh copy of the interpolating values array
 116  
      */
 117  
     public double[] getInterpolatingValues() {
 118  0
         double[] out = new double[y.length];
 119  0
         System.arraycopy(y, 0, out, 0, y.length);
 120  0
         return out;
 121  
     }
 122  
 
 123  
     /**
 124  
      * Returns a copy of the coefficients array.
 125  
      * <p>
 126  
      * Changes made to the returned copy will not affect the polynomial.
 127  
      * 
 128  
      * @return a fresh copy of the coefficients array
 129  
      */
 130  
     public double[] getCoefficients() {
 131  6
         if (!coefficientsComputed) {
 132  6
             computeCoefficients();
 133  
         }
 134  6
         double[] out = new double[coefficients.length];
 135  6
         System.arraycopy(coefficients, 0, out, 0, coefficients.length);
 136  6
         return out;
 137  
     }
 138  
 
 139  
     /**
 140  
      * Evaluate the Lagrange polynomial using 
 141  
      * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
 142  
      * Neville's Algorithm</a>. It takes O(N^2) time.
 143  
      * <p>
 144  
      * This function is made public static so that users can call it directly
 145  
      * without instantiating PolynomialFunctionLagrangeForm object.
 146  
      *
 147  
      * @param x the interpolating points array
 148  
      * @param y the interpolating values array
 149  
      * @param z the point at which the function value is to be computed
 150  
      * @return the function value
 151  
      * @throws FunctionEvaluationException if a runtime error occurs
 152  
      * @throws IllegalArgumentException if inputs are not valid
 153  
      */
 154  
     public static double evaluate(double x[], double y[], double z) throws
 155  
         FunctionEvaluationException, IllegalArgumentException {
 156  
 
 157  30
         int i, j, n, nearest = 0;
 158  
         double value, c[], d[], tc, td, divider, w, dist, min_dist;
 159  
 
 160  30
         verifyInterpolationArray(x, y);
 161  
 
 162  30
         n = x.length;
 163  30
         c = new double[n];
 164  30
         d = new double[n];
 165  30
         min_dist = Double.POSITIVE_INFINITY;
 166  158
         for (i = 0; i < n; i++) {
 167  
             // initialize the difference arrays
 168  128
             c[i] = y[i];
 169  128
             d[i] = y[i];
 170  
             // find out the abscissa closest to z
 171  128
             dist = Math.abs(z - x[i]);
 172  128
             if (dist < min_dist) {
 173  74
                 nearest = i;
 174  74
                 min_dist = dist;
 175  
             }
 176  
         }
 177  
 
 178  
         // initial approximation to the function value at z
 179  30
         value = y[nearest];
 180  
 
 181  122
         for (i = 1; i < n; i++) {
 182  330
             for (j = 0; j < n-i; j++) {
 183  238
                 tc = x[j] - z;
 184  238
                 td = x[i+j] - z;
 185  238
                 divider = x[j] - x[i+j];
 186  238
                 if (divider == 0.0) {
 187  
                     // This happens only when two abscissas are identical.
 188  2
                     throw new FunctionEvaluationException(z, 
 189  
                         "Identical abscissas cause division by zero: x[" +
 190  
                         i + "] = x[" + (i+j) + "] = " + x[i]);
 191  
                 }
 192  
                 // update the difference arrays
 193  236
                 w = (c[j+1] - d[j]) / divider;
 194  236
                 c[j] = tc * w;
 195  236
                 d[j] = td * w;
 196  
             }
 197  
             // sum up the difference terms to get the final value
 198  92
             if (nearest < 0.5*(n-i+1)) {
 199  34
                 value += c[nearest];    // fork down
 200  
             } else {
 201  58
                 nearest--;
 202  58
                 value += d[nearest];    // fork up
 203  
             }
 204  
         }
 205  
 
 206  28
         return value;
 207  
     }
 208  
 
 209  
     /**
 210  
      * Calculate the coefficients of Lagrange polynomial from the
 211  
      * interpolation data. It takes O(N^2) time.
 212  
      * <p>
 213  
      * Note this computation can be ill-conditioned. Use with caution
 214  
      * and only when it is necessary.
 215  
      *
 216  
      * @throws ArithmeticException if any abscissas coincide
 217  
      */
 218  
     protected void computeCoefficients() throws ArithmeticException {
 219  
         int i, j, n;
 220  
         double c[], tc[], d, t;
 221  
 
 222  6
         n = degree() + 1;
 223  6
         coefficients = new double[n];
 224  28
         for (i = 0; i < n; i++) {
 225  22
             coefficients[i] = 0.0;
 226  
         }
 227  
 
 228  
         // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
 229  6
         c = new double[n+1];
 230  6
         c[0] = 1.0;
 231  28
         for (i = 0; i < n; i++) {
 232  60
             for (j = i; j > 0; j--) {
 233  38
                 c[j] = c[j-1] - c[j] * x[i];
 234  
             }
 235  22
             c[0] *= (-x[i]);
 236  22
             c[i+1] = 1;
 237  
         }
 238  
 
 239  6
         tc = new double[n];
 240  28
         for (i = 0; i < n; i++) {
 241  
             // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
 242  22
             d = 1;
 243  120
             for (j = 0; j < n; j++) {
 244  98
                 if (i != j) {
 245  76
                     d *= (x[i] - x[j]);
 246  
                 }
 247  
             }
 248  22
             if (d == 0.0) {
 249  
                 // This happens only when two abscissas are identical.
 250  0
                 throw new ArithmeticException
 251  
                     ("Identical abscissas cause division by zero.");
 252  
             }
 253  22
             t = y[i] / d;
 254  
             // Lagrange polynomial is the sum of n terms, each of which is a
 255  
             // polynomial of degree n-1. tc[] are the coefficients of the i-th
 256  
             // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
 257  22
             tc[n-1] = c[n];     // actually c[n] = 1
 258  22
             coefficients[n-1] += t * tc[n-1];
 259  98
             for (j = n-2; j >= 0; j--) {
 260  76
                 tc[j] = c[j+1] + tc[j+1] * x[i];
 261  76
                 coefficients[j] += t * tc[j];
 262  
             }
 263  
         }
 264  
 
 265  6
         coefficientsComputed = true;
 266  6
     }
 267  
 
 268  
     /**
 269  
      * Verifies that the interpolation arrays are valid.
 270  
      * <p>
 271  
      * The interpolating points must be distinct. However it is not
 272  
      * verified here, it is checked in evaluate() and computeCoefficients().
 273  
      * 
 274  
      * @throws IllegalArgumentException if not valid
 275  
      * @see #evaluate(double[], double[], double)
 276  
      * @see #computeCoefficients()
 277  
      */
 278  
     protected static void verifyInterpolationArray(double x[], double y[]) throws
 279  
         IllegalArgumentException {
 280  
 
 281  58
         if (x.length < 2 || y.length < 2) {
 282  2
             throw new IllegalArgumentException
 283  
                 ("Interpolation requires at least two points.");
 284  
         }
 285  56
         if (x.length != y.length) {
 286  2
             throw new IllegalArgumentException
 287  
                 ("Abscissa and value arrays must have the same length.");
 288  
         }
 289  54
     }
 290  
 }