View Javadoc

    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 java.util.Arrays;
    20
    21  import org.apache.commons.math.FunctionEvaluationException;
    22
    23  /***
    24   * Represents a polynomial spline function.
    25   * <p>
    26   * A <strong>polynomial spline function</strong> consists of a set of
    27   * <i>interpolating polynomials</i> and an ascending array of domain
    28   * <i>knot points</i>, determining the intervals over which the spline function
    29   * is defined by the constituent polynomials.  The polynomials are assumed to
    30   * have been computed to match the values of another function at the knot
    31   * points.  The value consistency constraints are not currently enforced by
    32   * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
    33   * the polynomials and knot points passed to the constructor.
    34   * <p>
    35   * N.B.:  The polynomials in the <code>polynomials</code> property must be
    36   * centered on the knot points to compute the spline function values.  See below.
    37   * <p>
    38   * The domain of the polynomial spline function is
    39   * <code>[smallest knot, largest knot]</code>.  Attempts to evaluate the
    40   * function at values outside of this range generate IllegalArgumentExceptions.
    41   * <p>
    42   * The value of the polynomial spline function for an argument <code>x</code>
    43   * is computed as follows:
    44   * <ol>
    45   * <li>The knot array is searched to find the segment to which <code>x</code>
    46   * belongs.  If <code>x</code> is less than the smallest knot point or greater
    47   * than the largest one, an <code>IllegalArgumentException</code>
    48   * is thrown.</li>
    49   * <li> Let <code>j</code> be the index of the largest knot point that is less
    50   * than or equal to <code>x</code>.  The value returned is <br>
    51   * <code>polynomials[j](x - knot[j])</code></li></ol>
    52   *
    53   * @version $Revision: 348761 $ $Date: 2005-11-24 09:04:20 -0700 (Thu, 24 Nov 2005) $
    54   */
    55  public class PolynomialSplineFunction
    56      implements DifferentiableUnivariateRealFunction, Serializable {
    57
    58      /*** Serializable version identifier */
    59      private static final long serialVersionUID = 7011031166416885789L;
    60
    61      /*** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
    62      private double knots[];
    63
    64      /***
    65       * The polynomial functions that make up the spline.  The first element
    66       * determines the value of the spline over the first subinterval, the
    67       * second over the second, etc.   Spline function values are determined by
    68       * evaluating these functions at <code>(x - knot[i])</code> where i is the
    69       * knot segment to which x belongs.
    70       */
    71      private PolynomialFunction polynomials[] = null;
    72
    73      /***
    74       * Number of spline segments = number of polynomials
    75       *  = number of partition points - 1
    76       */
    77      private int n = 0;
    78
    79
    80      /***
    81       * Construct a polynomial spline function with the given segment delimiters
    82       * and interpolating polynomials.
    83       * <p>
    84       * The constructor copies both arrays and assigns the copies to the knots
    85       * and polynomials properties, respectively.
    86       *
    87       * @param knots spline segment interval delimiters
    88       * @param polynomials polynomial functions that make up the spline
    89       * @throws NullPointerException if either of the input arrays is null
    90       * @throws IllegalArgumentException if knots has length less than 2,
    91       * <code>polynomials.length != knots.length - 1 </code>, or the knots array
    92       * is not strictly increasing.
    93       *
    94       */
    95      public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
    96          if (knots.length < 2) {
    97              throw new IllegalArgumentException
    98                  ("Not enough knot values -- spline partition must have at least 2 points.");
    99          }
    100         if (knots.length - 1 != polynomials.length) {
    101             throw new IllegalArgumentException
    102             ("Number of polynomial interpolants must match the number of segments.");
    103         }
    104         if (!isStrictlyIncreasing(knots)) {
    105             throw new IllegalArgumentException
    106                 ("Knot values must be strictly increasing.");
    107         }
    108
    109         this.n = knots.length -1;
    110         this.knots = new double[n + 1];
    111         System.arraycopy(knots, 0, this.knots, 0, n + 1);
    112         this.polynomials = new PolynomialFunction[n];
    113         System.arraycopy(polynomials, 0, this.polynomials, 0, n);
    114     }
    115
    116     /***
    117      * Compute the value for the function.
    118      * <p>
    119      * Throws FunctionEvaluationException if v is outside of the domain of the
    120      * function.  The domain is [smallest knot, largest knot].
    121      * <p>
    122      * See {@link PolynomialSplineFunction} for details on the algorithm for
    123      * computing the value of the function.
    124      *
    125      * @param v the point for which the function value should be computed
    126      * @return the value
    127      * @throws FunctionEvaluationException if v is outside of the domain of
    128      * of the spline function (less than the smallest knot point or greater
    129      * than the largest knot point)
    130      */
    131     public double value(double v) throws FunctionEvaluationException {
    132         if (v < knots[0] || v > knots[n]) {
    133             throw new FunctionEvaluationException(v,"Argument outside domain");
    134         }
    135         int i = Arrays.binarySearch(knots, v);
    136         if (i < 0) {
    137             i = -i - 2;
    138         }
    139         //This will handle the case where v is the last knot value
    140         //There are only n-1 polynomials, so if v is the last knot
    141         //then we will use the last polynomial to calculate the value.
    142         if ( i >= polynomials.length ) {
    143             i--;
    144         }
    145         return polynomials[i].value(v - knots[i]);
    146     }
    147
    148     /***
    149      * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
    150      * @return  the derivative function
    151      */
    152     public UnivariateRealFunction derivative() {
    153         return polynomialSplineDerivative();
    154     }
    155
    156     /***
    157      * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
    158      *
    159      * @return  the derivative function
    160      */
    161     public PolynomialSplineFunction polynomialSplineDerivative() {
    162         PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
    163         for (int i = 0; i < n; i++) {
    164             derivativePolynomials[i] = polynomials[i].polynomialDerivative();
    165         }
    166         return new PolynomialSplineFunction(knots, derivativePolynomials);
    167     }
    168
    169     /***
    170      * Returns the number of spline segments = the number of polynomials
    171      * = the number of knot points - 1.
    172      *
    173      * @return the number of spline segments
    174      */
    175     public int getN() {
    176         return n;
    177     }
    178
    179     /***
    180      * Returns a copy of the interpolating polynomials array.
    181      * <p>
    182      * Returns a fresh copy of the array. Changes made to the copy will
    183      * not affect the polynomials property.
    184      *
    185      * @return the interpolating polynomials
    186      */
    187     public PolynomialFunction[] getPolynomials() {
    188         PolynomialFunction p[] = new PolynomialFunction[n];
    189         System.arraycopy(polynomials, 0, p, 0, n);
    190         return p;
    191     }
    192
    193     /***
    194      * Returns an array copy of the knot points.
    195      * <p>
    196      * Returns a fresh copy of the array. Changes made to the copy
    197      * will not affect the knots property.
    198      *
    199      * @return the knot points
    200      */
    201     public double[] getKnots() {
    202         double out[] = new double[n + 1];
    203         System.arraycopy(knots, 0, out, 0, n + 1);
    204         return out;
    205     }
    206
    207     /***
    208      * Determines if the given array is ordered in a strictly increasing
    209      * fashion.
    210      *
    211      * @param x the array to examine.
    212      * @return <code>true</code> if the elements in <code>x</code> are ordered
    213      * in a stricly increasing manner.  <code>false</code>, otherwise.
    214      */
    215     private static boolean isStrictlyIncreasing(double[] x) {
    216         for (int i = 1; i < x.length; ++i) {
    217             if (x[i - 1] >= x[i]) {
    218                 return false;
    219             }
    220         }
    221         return true;
    222     }
    223 }