You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by ps...@apache.org on 2004/04/02 23:16:24 UTC
cvs commit: jakarta-commons/math/src/java/org/apache/commons/math/analysis SplineInterpolator.java
psteitz 2004/04/02 13:16:24
Modified: math/src/java/org/apache/commons/math/analysis
SplineInterpolator.java
Log:
Fixed implementation, improved documentation.
Revision Changes Path
1.15 +88 -91 jakarta-commons/math/src/java/org/apache/commons/math/analysis/SplineInterpolator.java
Index: SplineInterpolator.java
===================================================================
RCS file: /home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/analysis/SplineInterpolator.java,v
retrieving revision 1.14
retrieving revision 1.15
diff -u -r1.14 -r1.15
--- SplineInterpolator.java 18 Feb 2004 03:24:19 -0000 1.14
+++ SplineInterpolator.java 2 Apr 2004 21:16:21 -0000 1.15
@@ -20,109 +20,106 @@
import java.io.Serializable;
/**
- * Computes a natural spline interpolation for the data set.
+ * Computes a natural (a.k.a. "free", "unclamped") cubic spline interpolation for the data set.
+ * <p>
+ * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction}
+ * consisting of n cubic polynomials, defined over the subintervals determined by the x values,
+ * x[0] < x[i] ... < x[n]. The x values are referred to as "knot points."
+ * <p>
+ * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest
+ * knot point and strictly less than the largest knot point is computed by finding the subinterval to which
+ * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where
+ * <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details.
+ * <p>
+ * The interpolating polynomials satisfy: <ol>
+ * <li>The value of the PolynomialSplineFunction at each of the input x values equals the
+ * corresponding y value.</li>
+ * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials
+ * "match up" at the knot points, as do their first and second derivatives).</li>
+ * </ol>
+ * <p>
+ * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires,
+ * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131.
*
* @version $Revision$ $Date$
*
*/
public class SplineInterpolator implements UnivariateRealInterpolator, Serializable {
- /** the natural spline coefficients. */
- private double[][] c = null;
-
+
/**
* Computes an interpolating function for the data set.
- * @param xval the arguments for the interpolation points
- * @param yval the values for the interpolation points
+ * @param x the arguments for the interpolation points
+ * @param y the values for the interpolation points
* @return a function which interpolates the data set
*/
- public UnivariateRealFunction interpolate(double[] xval, double[] yval) {
- if (xval.length != yval.length) {
+ public UnivariateRealFunction interpolate(double x[], double y[]) {
+ if (x.length != y.length) {
throw new IllegalArgumentException("Dataset arrays must have same length.");
}
-
- // TODO: What's this good for? Did I really write this???
- if (c == null) {
- // Number of intervals. The number of data points is N+1.
- int n = xval.length - 1;
- // Check whether the xval vector has ascending values.
- // TODO: Separation should be checked too (not implemented: which criteria?).
- for (int i = 0; i < n; i++) {
- if (xval[i] >= xval[i + 1]) {
- throw new IllegalArgumentException("Dataset must specify sorted, ascending x values.");
- }
- }
- // Vectors for the equation system. There are n-1 equations for the unknowns s[i] (1<=i<=N-1),
- // which are second order derivatives for the spline at xval[i]. At the end points, s[0]=s[N]=0.
- // Vector indices are offset by -1, except for the lower diagonal vector where the offset is -2. Layout of the equation system:
- // d[0]*s[1]+u[0]*s[2] = b[0]
- // l[0]*s[1]+d[1]*s[2]+u[1]*s[3] = b[1]
- // l[1]*s[2]+d[2]*s[3]+u[2]*s[4] = b[2]
- // ...
- // l[N-4]*s[N-3]+d[N-3]*s[N-2]+u[N-3]*s[N-1] = b[N-3]
- // l[N-3]*s[N-2]+d[N-2]*s[N-1] = b[N-2]
- // Vector b is the right hand side (RHS) of the system.
- double b[] = new double[n - 1];
- // Vector d is diagonal of the matrix and also holds the computed solution.
- double d[] = new double[n - 1];
- // Setup right hand side and diagonal.
- double dquot = (yval[1] - yval[0]) / (xval[1] - xval[0]);
- for (int i = 0; i < n - 1; i++) {
- double dquotNext =
- (yval[i + 2] - yval[i + 1]) / (xval[i + 2] - xval[i + 1]);
- b[i] = 6.0 * (dquotNext - dquot);
- d[i] = 2.0 * (xval[i + 2] - xval[i]);
- dquot = dquotNext;
- }
- // u[] and l[] (for the upper and lower diagonal respectively) are not
- // really needed, the computation is folded into the system solving loops.
- // Keep this for documentation purposes. The computation is folded into
- // the loops computing the solution.
- //double u[] = new double[n - 2]; // upper diagonal
- //double l[] = new double[n - 2]; // lower diagonal
- // Set up upper and lower diagonal. Keep the offsets in mind.
- //for (int i = 0; i < n - 2; i++) {
- // u[i] = xval[i + 2] - xval[i + 1];
- // l[i] = xval[i + 2] - xval[i + 1];
- //}
- // Solve the system: forward pass.
- for (int i = 0; i < n - 2; i++) {
- double delta = xval[i + 2] - xval[i + 1];
- double deltaquot = delta / d[i];
- d[i + 1] -= delta * deltaquot;
- b[i + 1] -= b[i] * deltaquot;
- }
- // Solve the system: backward pass.
- d[n - 2] = b[n - 2] / d[n - 2];
- for (int i = n - 3; i >= 0; i--) {
- d[i] = (b[i] - (xval[i + 2] - xval[i + 1]) * d[i + 1]) / d[i];
- }
- // Compute coefficients as usual polynomial coefficients.
- // Not the best with respect to roundoff on evaluation, but simple.
- c = new double[n][4];
- double delta = xval[1] - xval[0];
- c[0][3] = d[0] / delta / 6.0;
- c[0][2] = 0.0;
- c[0][1] = (yval[1] - yval[0]) / delta - d[0] * delta / 6.0;
- for (int i = 1; i < n - 2; i++) {
- delta = xval[i + 1] - xval[i];
- c[i][3] = (d[i] - d[i - 1]) / delta / 6.0;
- c[i][2] = d[i - 1] / 2.0;
- c[i][1] =
- (yval[i + 1] - yval[i]) / delta -
- (d[i] / 2.0 - d[i - 1]) * delta / 3.0;
- }
- delta = (xval[n] - xval[n - 1]);
- c[n - 1][3] = -d[n - 2] / delta / 6.0;
- c[n - 1][2] = d[n - 2] / 2.0;
- c[n - 1][1] =
- (yval[n] - yval[n - 1]) / delta - d[n - 2] * delta / 3.0;
- for (int i = 0; i < n; i++) {
- c[i][0] = yval[i];
+
+ if (x.length < 3) {
+ throw new IllegalArgumentException
+ ("At least 3 datapoints are required to compute a spline interpolant");
+ }
+
+ // Number of intervals. The number of data points is n + 1.
+ int n = x.length - 1;
+
+ for (int i = 0; i < n; i++) {
+ if (x[i] >= x[i + 1]) {
+ throw new IllegalArgumentException("Dataset x values must be strictly increasing.");
}
}
-
- // TODO: copy xval, unless copied in CubicSplineFunction constructor
- return new CubicSplineFunction(xval, c);
+
+ double h[] = new double[n];
+ for (int i = 0; i < n; i++) {
+ h[i] = x[i + 1] - x[i];
+ }
+
+ double alpha[] = new double[n];
+ for (int i = 1; i < n; i++) {
+ alpha[i] = 3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) /
+ (h[i - 1] * h[i]);
+ }
+
+ double l[] = new double[n + 1];
+ double mu[] = new double[n];
+ double z[] = new double[n + 1];
+ l[0] = 1d;
+ mu[0] = 0d;
+ z[0] = 0d;
+ for (int i = 1; i < n; i++) {
+ l[i] = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1];
+ mu[i] = h[i] / l[i];
+ z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
+ }
+
+ // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
+ double b[] = new double[n];
+ double c[] = new double[n + 1];
+ double d[] = new double[n];
+
+ l[n] = 1d;
+ z[n] = 0d;
+ c[n] = 0d;
+
+ for (int j = n -1; j >=0; j--) {
+ c[j] = z[j] - mu[j] * c[j + 1];
+ b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d;
+ d[j] = (c[j + 1] - c[j]) / (3d * h[j]);
+ }
+
+ PolynomialFunction polynomials[] = new PolynomialFunction[n];
+ double coefficients[] = new double[4];
+ for (int i = 0; i < n; i++) {
+ coefficients[0] = y[i];
+ coefficients[1] = b[i];
+ coefficients[2] = c[i];
+ coefficients[3] = d[i];
+ polynomials[i] = new PolynomialFunction(coefficients);
+ }
+
+ return new PolynomialSplineFunction(x, polynomials);
}
}
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org