Java API By Example, From Geeks To Geeks.

# Java > Open Source Codes > org > apache > commons > math > analysis > SplineInterpolator

 `1 /*2  * Copyright 2003-2004 The Apache Software Foundation.3  *4  * Licensed under the Apache License, Version 2.0 (the "License");5  * you may not use this file except in compliance with the License.6  * You may obtain a copy of the License at7  *8  * http://www.apache.org/licenses/LICENSE-2.09  *10  * Unless required by applicable law or agreed to in writing, software11  * distributed under the License is distributed on an "AS IS" BASIS,12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.13  * See the License for the specific language governing permissions and14  * limitations under the License.15  */16 package org.apache.commons.math.analysis;17 18 /**19  * Computes a natural (a.k.a. "free", "unclamped") cubic spline interpolation for the data set.20  *

21  * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction}22  * consisting of n cubic polynomials, defined over the subintervals determined by the x values, 23  * x[0] < x[i] ... < x[n]. The x values are referred to as "knot points."24  *

25  * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest26  * knot point and strictly less than the largest knot point is computed by finding the subinterval to which27  * x belongs and computing the value of the corresponding polynomial at x - x[i] where28  * i is the index of the subinterval. See {@link PolynomialSplineFunction} for more details.29  *

30  * The interpolating polynomials satisfy:

31  *
1. The value of the PolynomialSplineFunction at each of the input x values equals the 32  * corresponding y value.
2. 33  *
3. Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials 34  * "match up" at the knot points, as do their first and second derivatives).
4. 35  *
36  *

37  * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires, 38  * Numerical Analysis, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131.39  *40  * @version \$Revision\$ \$Date: 2005-02-26 05:11:52 -0800 (Sat, 26 Feb 2005) \$41  *42  */43 public class SplineInterpolator implements UnivariateRealInterpolator {44     45     /**46      * Computes an interpolating function for the data set.47      * @param x the arguments for the interpolation points48      * @param y the values for the interpolation points49      * @return a function which interpolates the data set50      */51     public UnivariateRealFunction interpolate(double x[], double y[]) {52         if (x.length != y.length) {53             throw new IllegalArgumentException ("Dataset arrays must have same length.");54         }55         56         if (x.length < 3) {57             throw new IllegalArgumentException 58                 ("At least 3 datapoints are required to compute a spline interpolant");59         }60         61         // Number of intervals. The number of data points is n + 1.62 int n = x.length - 1; 63         64         for (int i = 0; i < n; i++) {65             if (x[i] >= x[i + 1]) {66                 throw new IllegalArgumentException ("Dataset x values must be strictly increasing.");67             }68         }69         70         // Differences between knot points71 double h[] = new double[n];72         for (int i = 0; i < n; i++) {73             h[i] = x[i + 1] - x[i];74         }75         76         double mu[] = new double[n];77         double z[] = new double[n + 1];78         mu[0] = 0d;79         z[0] = 0d;80         double g = 0;81         for (int i = 1; i < n; i++) {82             g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1];83             mu[i] = h[i] / g;84             z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) /85                     (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g;86         }87        88         // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)89 double b[] = new double[n];90         double c[] = new double[n + 1];91         double d[] = new double[n];92         93         z[n] = 0d;94         c[n] = 0d;95         96         for (int j = n -1; j >=0; j--) {97             c[j] = z[j] - mu[j] * c[j + 1];98             b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d;99             d[j] = (c[j + 1] - c[j]) / (3d * h[j]);100         }101         102         PolynomialFunction polynomials[] = new PolynomialFunction[n];103         double coefficients[] = new double[4];104         for (int i = 0; i < n; i++) {105             coefficients[0] = y[i];106             coefficients[1] = b[i];107             coefficients[2] = c[i];108             coefficients[3] = d[i];109             polynomials[i] = new PolynomialFunction(coefficients);110         }111         112         return new PolynomialSplineFunction(x, polynomials);113     }114 115 }116 ` Popular Tags