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(
167 T[] x, T[] y)
168 throws MathIllegalArgumentException {
169
170 MathUtils.checkNotNull(x);
171 MathUtils.checkNotNull(y);
172 MathArrays.checkEqualLength(x, y);
173 if (x.length < 3) {
174 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
175 x.length, 3, true);
176 }
177
178 // Number of intervals. The number of data points is n + 1.
179 final int n = x.length - 1;
180
181 MathArrays.checkOrder(x);
182
183 // Differences between knot points
184 final Field<T> field = x[0].getField();
185 final T[] h = MathArrays.buildArray(field, n);
186 for (int i = 0; i < n; i++) {
187 h[i] = x[i + 1].subtract(x[i]);
188 }
189
190 final T[] mu = MathArrays.buildArray(field, n);
191 final T[] z = MathArrays.buildArray(field, n + 1);
192 mu[0] = field.getZero();
193 z[0] = field.getZero();
194 for (int i = 1; i < n; i++) {
195 final T g = x[i+1].subtract(x[i - 1]).multiply(2).subtract(h[i - 1].multiply(mu[i -1]));
196 mu[i] = h[i].divide(g);
197 z[i] = y[i + 1].multiply(h[i - 1]).
198 subtract(y[i].multiply(x[i + 1].subtract(x[i - 1]))).
199 add(y[i - 1].multiply(h[i])).
200 multiply(3).
201 divide(h[i - 1].multiply(h[i])).
202 subtract(h[i - 1].multiply(z[i - 1])).
203 divide(g);
204 }
205
206 // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
207 final T[] b = MathArrays.buildArray(field, n);
208 final T[] c = MathArrays.buildArray(field, n + 1);
209 final T[] d = MathArrays.buildArray(field, n);
210
211 z[n] = field.getZero();
212 c[n] = field.getZero();
213
214 for (int j = n -1; j >=0; j--) {
215 c[j] = z[j].subtract(mu[j].multiply(c[j + 1]));
216 b[j] = y[j + 1].subtract(y[j]).divide(h[j]).
217 subtract(h[j].multiply(c[j + 1].add(c[j]).add(c[j])).divide(3));
218 d[j] = c[j + 1].subtract(c[j]).divide(h[j].multiply(3));
219 }
220
221 @SuppressWarnings("unchecked")
222 final FieldPolynomialFunction<T>[] polynomials =
223 (FieldPolynomialFunction<T>[]) Array.newInstance(FieldPolynomialFunction.class, n);
224 final T[] coefficients = MathArrays.buildArray(field, 4);
225 for (int i = 0; i < n; i++) {
226 coefficients[0] = y[i];
227 coefficients[1] = b[i];
228 coefficients[2] = c[i];
229 coefficients[3] = d[i];
230 polynomials[i] = new FieldPolynomialFunction<>(coefficients);
231 }
232
233 return new FieldPolynomialSplineFunction<>(x, polynomials);
234 }
235
236 }