Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 1 | /* |
| 2 | * Licensed to the Apache Software Foundation (ASF) under one or more |
| 3 | * contributor license agreements. See the NOTICE file distributed with |
| 4 | * this work for additional information regarding copyright ownership. |
| 5 | * The ASF licenses this file to You under the Apache License, Version 2.0 |
| 6 | * (the "License"); you may not use this file except in compliance with |
| 7 | * the License. You may obtain a copy of the License at |
| 8 | * |
| 9 | * http://www.apache.org/licenses/LICENSE-2.0 |
| 10 | * |
| 11 | * Unless required by applicable law or agreed to in writing, software |
| 12 | * distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | * See the License for the specific language governing permissions and |
| 15 | * limitations under the License. |
| 16 | */ |
| 17 | package org.apache.commons.math.analysis.interpolation; |
| 18 | |
| 19 | import org.apache.commons.math.exception.DimensionMismatchException; |
| 20 | import org.apache.commons.math.exception.util.LocalizedFormats; |
| 21 | import org.apache.commons.math.exception.NumberIsTooSmallException; |
| 22 | import org.apache.commons.math.analysis.polynomials.PolynomialFunction; |
| 23 | import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; |
| 24 | import org.apache.commons.math.util.MathUtils; |
| 25 | |
| 26 | /** |
| 27 | * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set. |
| 28 | * <p> |
| 29 | * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction} |
| 30 | * consisting of n cubic polynomials, defined over the subintervals determined by the x values, |
| 31 | * x[0] < x[i] ... < x[n]. The x values are referred to as "knot points."</p> |
| 32 | * <p> |
| 33 | * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest |
| 34 | * knot point and strictly less than the largest knot point is computed by finding the subinterval to which |
| 35 | * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where |
| 36 | * <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details. |
| 37 | * </p> |
| 38 | * <p> |
| 39 | * The interpolating polynomials satisfy: <ol> |
| 40 | * <li>The value of the PolynomialSplineFunction at each of the input x values equals the |
| 41 | * corresponding y value.</li> |
| 42 | * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials |
| 43 | * "match up" at the knot points, as do their first and second derivatives).</li> |
| 44 | * </ol></p> |
| 45 | * <p> |
| 46 | * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires, |
| 47 | * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131. |
| 48 | * </p> |
| 49 | * |
| 50 | * @version $Revision: 983921 $ $Date: 2010-08-10 12:46:06 +0200 (mar. 10 août 2010) $ |
| 51 | * |
| 52 | */ |
| 53 | public class SplineInterpolator implements UnivariateRealInterpolator { |
| 54 | |
| 55 | /** |
| 56 | * Computes an interpolating function for the data set. |
| 57 | * @param x the arguments for the interpolation points |
| 58 | * @param y the values for the interpolation points |
| 59 | * @return a function which interpolates the data set |
| 60 | * @throws DimensionMismatchException if {@code x} and {@code y} |
| 61 | * have different sizes. |
| 62 | * @throws org.apache.commons.math.exception.NonMonotonousSequenceException |
| 63 | * if {@code x} is not sorted in strict increasing order. |
| 64 | * @throws NumberIsTooSmallException if the size of {@code x} is smaller |
| 65 | * than 3. |
| 66 | */ |
| 67 | public PolynomialSplineFunction interpolate(double x[], double y[]) { |
| 68 | if (x.length != y.length) { |
| 69 | throw new DimensionMismatchException(x.length, y.length); |
| 70 | } |
| 71 | |
| 72 | if (x.length < 3) { |
| 73 | throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, |
| 74 | x.length, 3, true); |
| 75 | } |
| 76 | |
| 77 | // Number of intervals. The number of data points is n + 1. |
| 78 | int n = x.length - 1; |
| 79 | |
| 80 | MathUtils.checkOrder(x); |
| 81 | |
| 82 | // Differences between knot points |
| 83 | double h[] = new double[n]; |
| 84 | for (int i = 0; i < n; i++) { |
| 85 | h[i] = x[i + 1] - x[i]; |
| 86 | } |
| 87 | |
| 88 | double mu[] = new double[n]; |
| 89 | double z[] = new double[n + 1]; |
| 90 | mu[0] = 0d; |
| 91 | z[0] = 0d; |
| 92 | double g = 0; |
| 93 | for (int i = 1; i < n; i++) { |
| 94 | g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1]; |
| 95 | mu[i] = h[i] / g; |
| 96 | z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) / |
| 97 | (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g; |
| 98 | } |
| 99 | |
| 100 | // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) |
| 101 | double b[] = new double[n]; |
| 102 | double c[] = new double[n + 1]; |
| 103 | double d[] = new double[n]; |
| 104 | |
| 105 | z[n] = 0d; |
| 106 | c[n] = 0d; |
| 107 | |
| 108 | for (int j = n -1; j >=0; j--) { |
| 109 | c[j] = z[j] - mu[j] * c[j + 1]; |
| 110 | b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d; |
| 111 | d[j] = (c[j + 1] - c[j]) / (3d * h[j]); |
| 112 | } |
| 113 | |
| 114 | PolynomialFunction polynomials[] = new PolynomialFunction[n]; |
| 115 | double coefficients[] = new double[4]; |
| 116 | for (int i = 0; i < n; i++) { |
| 117 | coefficients[0] = y[i]; |
| 118 | coefficients[1] = b[i]; |
| 119 | coefficients[2] = c[i]; |
| 120 | coefficients[3] = d[i]; |
| 121 | polynomials[i] = new PolynomialFunction(coefficients); |
| 122 | } |
| 123 | |
| 124 | return new PolynomialSplineFunction(x, polynomials); |
| 125 | } |
| 126 | |
| 127 | } |