1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
140
141
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 }