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 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 }