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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.interpolation;
23
24 import java.lang.reflect.Array;
25
26 import org.hipparchus.Field;
27 import org.hipparchus.CalculusFieldElement;
28 import org.hipparchus.analysis.polynomials.FieldPolynomialFunction;
29 import org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction;
30 import org.hipparchus.analysis.polynomials.PolynomialFunction;
31 import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
32 import org.hipparchus.exception.LocalizedCoreFormats;
33 import org.hipparchus.exception.MathIllegalArgumentException;
34 import org.hipparchus.util.MathArrays;
35 import org.hipparchus.util.MathUtils;
36
37 /**
38 * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set.
39 * <p>
40 * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction}
41 * consisting of n cubic polynomials, defined over the subintervals determined by the x values,
42 * {@code x[0] < x[i] ... < x[n].} The x values are referred to as "knot points."</p>
43 * <p>
44 * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest
45 * knot point and strictly less than the largest knot point is computed by finding the subinterval to which
46 * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where
47 * <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details.
48 * </p>
49 * <p>
50 * The interpolating polynomials satisfy:
51 * </p>
52 * <ol>
53 * <li>The value of the PolynomialSplineFunction at each of the input x values equals the
54 * corresponding y value.</li>
55 * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials
56 * "match up" at the knot points, as do their first and second derivatives).</li>
57 * </ol>
58 * <p>
59 * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires,
60 * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131.
61 * </p>
62 *
63 */
64 public class SplineInterpolator implements UnivariateInterpolator, FieldUnivariateInterpolator {
65
66 /** Empty constructor.
67 * <p>
68 * This constructor is not strictly necessary, but it prevents spurious
69 * javadoc warnings with JDK 18 and later.
70 * </p>
71 * @since 3.0
72 */
73 public SplineInterpolator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
74 // nothing to do
75 }
76
77 /**
78 * Computes an interpolating function for the data set.
79 * @param x the arguments for the interpolation points
80 * @param y the values for the interpolation points
81 * @return a function which interpolates the data set
82 * @throws MathIllegalArgumentException if {@code x} and {@code y}
83 * have different sizes.
84 * @throws MathIllegalArgumentException if {@code x} is not sorted in
85 * strict increasing order.
86 * @throws MathIllegalArgumentException if the size of {@code x} is smaller
87 * than 3.
88 */
89 @Override
90 public PolynomialSplineFunction interpolate(double x[], double y[])
91 throws MathIllegalArgumentException {
92
93 MathUtils.checkNotNull(x);
94 MathUtils.checkNotNull(y);
95 MathArrays.checkEqualLength(x, y);
96 if (x.length < 3) {
97 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
98 x.length, 3, true);
99 }
100
101 // Number of intervals. The number of data points is n + 1.
102 final int n = x.length - 1;
103
104 MathArrays.checkOrder(x);
105
106 // Differences between knot points
107 final double h[] = new double[n];
108 for (int i = 0; i < n; i++) {
109 h[i] = x[i + 1] - x[i];
110 }
111
112 final double mu[] = new double[n];
113 final double z[] = new double[n + 1];
114 mu[0] = 0d;
115 z[0] = 0d;
116 double g;
117 for (int i = 1; i < n; i++) {
118 g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1];
119 mu[i] = h[i] / g;
120 z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) /
121 (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g;
122 }
123
124 // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
125 final double b[] = new double[n];
126 final double c[] = new double[n + 1];
127 final double d[] = new double[n];
128
129 z[n] = 0d;
130 c[n] = 0d;
131
132 for (int j = n -1; j >=0; j--) {
133 c[j] = z[j] - mu[j] * c[j + 1];
134 b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d;
135 d[j] = (c[j + 1] - c[j]) / (3d * h[j]);
136 }
137
138 final PolynomialFunction polynomials[] = new PolynomialFunction[n];
139 final double coefficients[] = new double[4];
140 for (int i = 0; i < n; i++) {
141 coefficients[0] = y[i];
142 coefficients[1] = b[i];
143 coefficients[2] = c[i];
144 coefficients[3] = d[i];
145 polynomials[i] = new PolynomialFunction(coefficients);
146 }
147
148 return new PolynomialSplineFunction(x, polynomials);
149 }
150
151 /**
152 * Computes an interpolating function for the data set.
153 * @param x the arguments for the interpolation points
154 * @param y the values for the interpolation points
155 * @param <T> the type of the field elements
156 * @return a function which interpolates the data set
157 * @throws MathIllegalArgumentException if {@code x} and {@code y}
158 * have different sizes.
159 * @throws MathIllegalArgumentException if {@code x} is not sorted in
160 * strict increasing order.
161 * @throws MathIllegalArgumentException if the size of {@code x} is smaller
162 * than 3.
163 * @since 1.5
164 */
165 @Override
166 public <T extends CalculusFieldElement<T>> FieldPolynomialSplineFunction<T> interpolate(T x[], T y[])
167 throws MathIllegalArgumentException {
168
169 MathUtils.checkNotNull(x);
170 MathUtils.checkNotNull(y);
171 MathArrays.checkEqualLength(x, y);
172 if (x.length < 3) {
173 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
174 x.length, 3, true);
175 }
176
177 // Number of intervals. The number of data points is n + 1.
178 final int n = x.length - 1;
179
180 MathArrays.checkOrder(x);
181
182 // Differences between knot points
183 final Field<T> field = x[0].getField();
184 final T h[] = MathArrays.buildArray(field, n);
185 for (int i = 0; i < n; i++) {
186 h[i] = x[i + 1].subtract(x[i]);
187 }
188
189 final T mu[] = MathArrays.buildArray(field, n);
190 final T z[] = MathArrays.buildArray(field, n + 1);
191 mu[0] = field.getZero();
192 z[0] = field.getZero();
193 for (int i = 1; i < n; i++) {
194 final T g = x[i+1].subtract(x[i - 1]).multiply(2).subtract(h[i - 1].multiply(mu[i -1]));
195 mu[i] = h[i].divide(g);
196 z[i] = y[i + 1].multiply(h[i - 1]).
197 subtract(y[i].multiply(x[i + 1].subtract(x[i - 1]))).
198 add(y[i - 1].multiply(h[i])).
199 multiply(3).
200 divide(h[i - 1].multiply(h[i])).
201 subtract(h[i - 1].multiply(z[i - 1])).
202 divide(g);
203 }
204
205 // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
206 final T b[] = MathArrays.buildArray(field, n);
207 final T c[] = MathArrays.buildArray(field, n + 1);
208 final T d[] = MathArrays.buildArray(field, n);
209
210 z[n] = field.getZero();
211 c[n] = field.getZero();
212
213 for (int j = n -1; j >=0; j--) {
214 c[j] = z[j].subtract(mu[j].multiply(c[j + 1]));
215 b[j] = y[j + 1].subtract(y[j]).divide(h[j]).
216 subtract(h[j].multiply(c[j + 1].add(c[j]).add(c[j])).divide(3));
217 d[j] = c[j + 1].subtract(c[j]).divide(h[j].multiply(3));
218 }
219
220 @SuppressWarnings("unchecked")
221 final FieldPolynomialFunction<T> polynomials[] =
222 (FieldPolynomialFunction<T>[]) Array.newInstance(FieldPolynomialFunction.class, n);
223 final T coefficients[] = MathArrays.buildArray(field, 4);
224 for (int i = 0; i < n; i++) {
225 coefficients[0] = y[i];
226 coefficients[1] = b[i];
227 coefficients[2] = c[i];
228 coefficients[3] = d[i];
229 polynomials[i] = new FieldPolynomialFunction<>(coefficients);
230 }
231
232 return new FieldPolynomialSplineFunction<>(x, polynomials);
233 }
234
235 }