View Javadoc

    1   /*
    2    * Copyright 2003-2004 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 or equal to 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: 1.9 $ $Date: 2004/07/22 02:34:25 $
    54   */
    55  public class PolynomialSplineFunction implements UnivariateRealFunction, Serializable {
    56
    57      /*** Serializable version identifier */
    58      static final long serialVersionUID = 7011031166416885789L;
    59
    60      /*** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
    61      private double knots[];
    62
    63      /***
    64       * The polynomial functions that make up the spline.  The first element
    65       * determines the value of the spline over the first subinterval, the
    66       * second over the second, etc.   Spline function values are determined by
    67       * evaluating these functions at <code>(x - knot[i])</code> where i is the
    68       * knot segment to which x belongs.
    69       */
    70      private PolynomialFunction polynomials[] = null;
    71
    72      /***
    73       * Number of spline segments = number of polynomials
    74       *  = number of partition points - 1
    75       */
    76      private int n = 0;
    77
    78
    79      /***
    80       * Construct a polynomial spline function with the given segment delimiters
    81       * and interpolating polynomials.
    82       * <p>
    83       * The constructor copies both arrays and assigns the copies to the knots
    84       * and polynomials properties, respectively.
    85       *
    86       * @param knots spline segment interval delimiters
    87       * @param polynomials polynomial functions that make up the spline
    88       * @throws NullPointerException if either of the input arrays is null
    89       * @throws IllegalArgumentException if knots has length less than 2,
    90       * <code>polynomials.length != knots.length - 1 </code>, or the knots array
    91       * is not strictly increasing.
    92       *
    93       */
    94      public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
    95          if (knots.length < 2) {
    96              throw new IllegalArgumentException
    97              	("Not enough knot values -- spline partition must have at least 2 points.");
    98          }
    99          if (knots.length - 1 != polynomials.length) {
    100             throw new IllegalArgumentException
    101             ("Number of polynomial interpolants must match the number of segments.");
    102         }
    103         if (!isStrictlyIncreasing(knots)) {
    104             throw new IllegalArgumentException
    105                 ("Knot values must be strictly increasing.");
    106         }
    107
    108         this.n = knots.length -1;
    109         this.knots = new double[n + 1];
    110         System.arraycopy(knots, 0, this.knots, 0, n + 1);
    111         this.polynomials = new PolynomialFunction[n];
    112         System.arraycopy(polynomials, 0, this.polynomials, 0, n);
    113     }
    114
    115     /***
    116      * Compute the value for the function.
    117      * <p>
    118      * Throws FunctionEvaluationException if v is outside of the domain of the
    119      * function.  The domain is [smallest knot, largest knot).
    120      * <p>
    121      * See {@link PolynomialSplineFunction} for details on the algorithm for
    122      * computing the value of the function.
    123      *
    124      * @param v the point for which the function value should be computed
    125      * @return the value
    126      * @throws FunctionEvaluationException if v is outside of the domain of
    127      * of the spline function (less than the smallest knot point or greater
    128      * than or equal to the largest knot point)
    129      */
    130     public double value(double v) throws FunctionEvaluationException {
    131         if (v < knots[0] || v >= knots[n]) {
    132             throw new FunctionEvaluationException(v,"Argument outside domain");
    133         }
    134         int i = Arrays.binarySearch(knots, v);
    135         if (i < 0) {
    136             i = -i - 2;
    137         }
    138         return polynomials[i].value(v - knots[i]);
    139     }
    140
    141     /***
    142      * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
    143      * @return  the derivative function
    144      */
    145     public UnivariateRealFunction derivative() {
    146         return polynomialSplineDerivative();
    147     }
    148
    149     /***
    150      * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
    151      *
    152      * @return  the derivative function
    153      */
    154     public PolynomialSplineFunction polynomialSplineDerivative() {
    155         PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
    156         for (int i = 0; i < n; i++) {
    157             derivativePolynomials[i] = polynomials[i].polynomialDerivative();
    158         }
    159         return new PolynomialSplineFunction(knots, derivativePolynomials);
    160     }
    161
    162     /***
    163      * Returns the number of spline segments = the number of polynomials
    164      * = the number of knot points - 1.
    165      *
    166      * @return the number of spline segments
    167      */
    168     public int getN() {
    169         return n;
    170     }
    171
    172     /***
    173      * Returns a copy of the interpolating polynomials array.
    174      * <p>
    175      * Returns a fresh copy of the array. Changes made to the copy will
    176      * not affect the polynomials property.
    177      *
    178      * @return the interpolating polynomials
    179      */
    180     public PolynomialFunction[] getPolynomials() {
    181         PolynomialFunction p[] = new PolynomialFunction[n];
    182         System.arraycopy(polynomials, 0, p, 0, n);
    183         return p;
    184     }
    185
    186     /***
    187      * Returns an array copy of the knot points.
    188      * <p>
    189      * Returns a fresh copy of the array. Changes made to the copy
    190      * will not affect the knots property.
    191      *
    192      * @return the knot points
    193      */
    194     public double[] getKnots() {
    195         double out[] = new double[n + 1];
    196         System.arraycopy(knots, 0, out, 0, n + 1);
    197         return out;
    198     }
    199
    200     /***
    201      * Determines if the given array is ordered in a strictly increasing
    202      * fashion.
    203      *
    204      * @param x the array to examine.
    205      * @return <code>true</code> if the elements in <code>x</code> are ordered
    206      * in a stricly increasing manner.  <code>false</code>, otherwise.
    207      */
    208     private static boolean isStrictlyIncreasing(double[] x) {
    209         for (int i = 1; i < x.length; ++i) {
    210             if (x[i - 1] >= x[i]) {
    211                 return false;
    212             }
    213         }
    214         return true;
    215     }
    216 }