001: /*
002: * Copyright 2003-2005 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.analysis;
017:
018: import java.io.Serializable;
019: import java.util.Arrays;
020:
021: import org.apache.commons.math.FunctionEvaluationException;
022:
023: /**
024: * Represents a polynomial spline function.
025: * <p>
026: * A <strong>polynomial spline function</strong> consists of a set of
027: * <i>interpolating polynomials</i> and an ascending array of domain
028: * <i>knot points</i>, determining the intervals over which the spline function
029: * is defined by the constituent polynomials. The polynomials are assumed to
030: * have been computed to match the values of another function at the knot
031: * points. The value consistency constraints are not currently enforced by
032: * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
033: * the polynomials and knot points passed to the constructor.
034: * <p>
035: * N.B.: The polynomials in the <code>polynomials</code> property must be
036: * centered on the knot points to compute the spline function values. See below.
037: * <p>
038: * The domain of the polynomial spline function is
039: * <code>[smallest knot, largest knot]</code>. Attempts to evaluate the
040: * function at values outside of this range generate IllegalArgumentExceptions.
041: * <p>
042: * The value of the polynomial spline function for an argument <code>x</code>
043: * is computed as follows:
044: * <ol>
045: * <li>The knot array is searched to find the segment to which <code>x</code>
046: * belongs. If <code>x</code> is less than the smallest knot point or greater
047: * than the largest one, an <code>IllegalArgumentException</code>
048: * is thrown.</li>
049: * <li> Let <code>j</code> be the index of the largest knot point that is less
050: * than or equal to <code>x</code>. The value returned is <br>
051: * <code>polynomials[j](x - knot[j])</code></li></ol>
052: *
053: * @version $Revision: 348761 $ $Date: 2005-11-24 09:04:20 -0700 (Thu, 24 Nov 2005) $
054: */
055: public class PolynomialSplineFunction implements
056: DifferentiableUnivariateRealFunction, Serializable {
057:
058: /** Serializable version identifier */
059: private static final long serialVersionUID = 7011031166416885789L;
060:
061: /** Spline segment interval delimiters (knots). Size is n+1 for n segments. */
062: private double knots[];
063:
064: /**
065: * The polynomial functions that make up the spline. The first element
066: * determines the value of the spline over the first subinterval, the
067: * second over the second, etc. Spline function values are determined by
068: * evaluating these functions at <code>(x - knot[i])</code> where i is the
069: * knot segment to which x belongs.
070: */
071: private PolynomialFunction polynomials[] = null;
072:
073: /**
074: * Number of spline segments = number of polynomials
075: * = number of partition points - 1
076: */
077: private int n = 0;
078:
079: /**
080: * Construct a polynomial spline function with the given segment delimiters
081: * and interpolating polynomials.
082: * <p>
083: * The constructor copies both arrays and assigns the copies to the knots
084: * and polynomials properties, respectively.
085: *
086: * @param knots spline segment interval delimiters
087: * @param polynomials polynomial functions that make up the spline
088: * @throws NullPointerException if either of the input arrays is null
089: * @throws IllegalArgumentException if knots has length less than 2,
090: * <code>polynomials.length != knots.length - 1 </code>, or the knots array
091: * is not strictly increasing.
092: *
093: */
094: public PolynomialSplineFunction(double knots[],
095: PolynomialFunction polynomials[]) {
096: if (knots.length < 2) {
097: throw new IllegalArgumentException(
098: "Not enough knot values -- spline partition must have at least 2 points.");
099: }
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,
134: "Argument outside domain");
135: }
136: int i = Arrays.binarySearch(knots, v);
137: if (i < 0) {
138: i = -i - 2;
139: }
140: //This will handle the case where v is the last knot value
141: //There are only n-1 polynomials, so if v is the last knot
142: //then we will use the last polynomial to calculate the value.
143: if (i >= polynomials.length) {
144: i--;
145: }
146: return polynomials[i].value(v - knots[i]);
147: }
148:
149: /**
150: * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
151: * @return the derivative function
152: */
153: public UnivariateRealFunction derivative() {
154: return polynomialSplineDerivative();
155: }
156:
157: /**
158: * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
159: *
160: * @return the derivative function
161: */
162: public PolynomialSplineFunction polynomialSplineDerivative() {
163: PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
164: for (int i = 0; i < n; i++) {
165: derivativePolynomials[i] = polynomials[i]
166: .polynomialDerivative();
167: }
168: return new PolynomialSplineFunction(knots,
169: derivativePolynomials);
170: }
171:
172: /**
173: * Returns the number of spline segments = the number of polynomials
174: * = the number of knot points - 1.
175: *
176: * @return the number of spline segments
177: */
178: public int getN() {
179: return n;
180: }
181:
182: /**
183: * Returns a copy of the interpolating polynomials array.
184: * <p>
185: * Returns a fresh copy of the array. Changes made to the copy will
186: * not affect the polynomials property.
187: *
188: * @return the interpolating polynomials
189: */
190: public PolynomialFunction[] getPolynomials() {
191: PolynomialFunction p[] = new PolynomialFunction[n];
192: System.arraycopy(polynomials, 0, p, 0, n);
193: return p;
194: }
195:
196: /**
197: * Returns an array copy of the knot points.
198: * <p>
199: * Returns a fresh copy of the array. Changes made to the copy
200: * will not affect the knots property.
201: *
202: * @return the knot points
203: */
204: public double[] getKnots() {
205: double out[] = new double[n + 1];
206: System.arraycopy(knots, 0, out, 0, n + 1);
207: return out;
208: }
209:
210: /**
211: * Determines if the given array is ordered in a strictly increasing
212: * fashion.
213: *
214: * @param x the array to examine.
215: * @return <code>true</code> if the elements in <code>x</code> are ordered
216: * in a stricly increasing manner. <code>false</code>, otherwise.
217: */
218: private static boolean isStrictlyIncreasing(double[] x) {
219: for (int i = 1; i < x.length; ++i) {
220: if (x[i - 1] >= x[i]) {
221: return false;
222: }
223: }
224: return true;
225: }
226: }
|