commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject cvs commit: jakarta-commons/math/src/java/org/apache/commons/math/analysis SplineInterpolator.java
Date Fri, 02 Apr 2004 21:16:24 GMT
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


Mime
View raw message