1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
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 at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * 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 and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.interpolation;
17
18 /**
19 * Computes cubic spline interpolation.
20 * Accuracy of this interpolator worsens as the polynomial degree to interpolate increases.
21 * Typically, up to degree 3 results are reasonable.
22 */
23 public class CubicSplineInterpolator extends BaseInterpolator {
24
25 /**
26 * Length of x's and y's to take into account.
27 */
28 private static final int M = 2;
29
30 private static final double YP1 = 1.e99;
31
32 private static final double YPN = 1.e99;
33
34 private final double[] y2;
35
36 /**
37 * Constructor with x and y vectors, and values of the first derivative at the endpoints.
38 *
39 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
40 * decreasing)
41 * @param y y values to interpolate to.
42 * @param yp1 1st derivative at lowest endpoint.
43 * @param ypn 1st derivative at highest endpoint
44 */
45 public CubicSplineInterpolator(final double[] x, final double[] y, final double yp1, final double ypn) {
46 super(x, y, M);
47 y2 = new double[x.length];
48 setY2(x, y, yp1, ypn);
49 }
50
51 /**
52 * Constructor with x and y vectors and default values for first derivative at the endpoints.
53 *
54 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
55 * decreasing)
56 * @param y y values to interpolate to.
57 */
58 public CubicSplineInterpolator(final double[] x, final double[] y) {
59 this(x, y, YP1, YPN);
60 }
61
62 /**
63 * Actual interpolation method.
64 *
65 * @param jl index where value x to be interpolated in located in the array of xx.
66 * @param x value to obtain interpolation for.
67 * @return interpolated value.
68 * @throws InterpolationException if interpolation fails.
69 */
70 @Override
71 public double rawinterp(int jl, double x) throws InterpolationException {
72 // Given a value x, and using pointers to data xx and yy, and the stored vector of second
73 // derivatives y2, this routine this cubic spline interpolated value y
74 final var khi = jl + 1;
75
76 final var h = xx[khi] - xx[jl];
77 if (h == 0.0) {
78 // The xa's must be distinct
79 throw new InterpolationException();
80 }
81
82 // Cubic spline polynomial is now evaluated
83 final var a = (xx[khi] - x) / h;
84 final var b = (x - xx[jl]) / h;
85 return a * yy[jl] + b * yy[khi] + ((a * a * a - a) * y2[jl] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
86 }
87
88 /**
89 * This method stores an array y2[0..n-1] with second derivatives of the interpolating function
90 * at the tabulated points pointed to by xv, using function values pointed to by yv. If yp1
91 * and/or ypn are equal to 1e99 or larger, the routine is signaled to set the corresponding
92 * boundary condition for a natural spline, with zero second derivative on that boundary;
93 * otherwise, they are the values of the first derivatives at the endpoints.
94 *
95 * @param xv x values to interpolate to. Values in x must be monotonic (either increasing or
96 * decreasing)
97 * @param yv y values to interpolate to.
98 * @param yp1 1st derivative at lowest endpoint.
99 * @param ypn 1st derivative at highest endpoint
100 */
101 private void setY2(final double[] xv, final double[] yv, final double yp1, final double ypn) {
102 int i;
103 int k;
104 double p;
105 final double qn;
106 double sig;
107 final double un;
108 final var n = y2.length;
109 final var u = new double[n - 1];
110
111 if (yp1 > 0.99e99) {
112 // The lower boundary condition is set either to be "natural"
113 y2[0] = u[0] = 0.0;
114 } else {
115 // or else to have a specified first derivative
116 y2[0] = -0.5;
117 u[0] = (3.0 / (xv[1] - xv[0])) * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);
118 }
119
120 for (i = 1; i < n - 1; i++) {
121 // This is the decomposition loop of the tri-diagonal algorithm. y2 and "u" are used for
122 // temporary storage of the decomposed factors
123 sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
124 p = sig * y2[i - 1] + 2.0;
125 y2[i] = (sig - 1.0) / p;
126 u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]) - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
127 u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
128 }
129
130 if (ypn > 0.99e99) {
131 // The upper boundary condition is set either to be "natural"
132 qn = un = 0.0;
133 } else {
134 // or else to have a specified first derivative
135 qn = 0.5;
136 un = (3.0 / (xv[n - 1] - xv[n - 2])) * (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));
137 }
138 y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
139 for (k = n - 2; k >= 0; k--) {
140 // This is the back-substitution loop of the tri-diagonal algorithm
141 y2[k] = y2[k] * y2[k + 1] + u[k];
142 }
143 }
144 }