Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed implementation, improved documentation. 

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@141145 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz 2004-04-02 21:16:21 +00:00
parent 0a46550083
commit 36fcb5c54d
1 changed files with 88 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

179
src/java/org/apache/commons/math/analysis/SplineInterpolator.java
View File

@ -20,109 +20,106 @@ package org.apache.commons.math.analysis;
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: 1.14 $ $Date: 2004/02/18 03:24:19 $
* @version $Revision: 1.15 $ $Date: 2004/04/02 21:16:21 $
*
*/
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);
}
}