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.math4.legacy.analysis.interpolation;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.apache.commons.math4.legacy.core.FieldElement;
23 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
24 import org.apache.commons.math4.legacy.exception.MathArithmeticException;
25 import org.apache.commons.math4.legacy.exception.NoDataException;
26 import org.apache.commons.math4.legacy.exception.NullArgumentException;
27 import org.apache.commons.math4.legacy.exception.ZeroException;
28 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
29 import org.apache.commons.math4.legacy.core.MathArrays;
30
31 /** Polynomial interpolator using both sample values and sample derivatives.
32 * <p>
33 * The interpolation polynomials match all sample points, including both values
34 * and provided derivatives. There is one polynomial for each component of
35 * the values vector. All polynomials have the same degree. The degree of the
36 * polynomials depends on the number of points and number of derivatives at each
37 * point. For example the interpolation polynomials for n sample points without
38 * any derivatives all have degree n-1. The interpolation polynomials for n
39 * sample points with the two extreme points having value and first derivative
40 * and the remaining points having value only all have degree n+1. The
41 * interpolation polynomial for n sample points with value, first and second
42 * derivative for all points all have degree 3n-1.
43 * </p>
44 *
45 * @param <T> Type of the field elements.
46 *
47 * @since 3.2
48 */
49 public class FieldHermiteInterpolator<T extends FieldElement<T>> {
50
51 /** Sample abscissae. */
52 private final List<T> abscissae;
53
54 /** Top diagonal of the divided differences array. */
55 private final List<T[]> topDiagonal;
56
57 /** Bottom diagonal of the divided differences array. */
58 private final List<T[]> bottomDiagonal;
59
60 /** Create an empty interpolator.
61 */
62 public FieldHermiteInterpolator() {
63 this.abscissae = new ArrayList<>();
64 this.topDiagonal = new ArrayList<>();
65 this.bottomDiagonal = new ArrayList<>();
66 }
67
68 /** Add a sample point.
69 * <p>
70 * This method must be called once for each sample point. It is allowed to
71 * mix some calls with values only with calls with values and first
72 * derivatives.
73 * </p>
74 * <p>
75 * The point abscissae for all calls <em>must</em> be different.
76 * </p>
77 * @param x abscissa of the sample point
78 * @param value value and derivatives of the sample point
79 * (if only one row is passed, it is the value, if two rows are
80 * passed the first one is the value and the second the derivative
81 * and so on)
82 * @exception ZeroException if the abscissa difference between added point
83 * and a previous point is zero (i.e. the two points are at same abscissa)
84 * @exception MathArithmeticException if the number of derivatives is larger
85 * than 20, which prevents computation of a factorial
86 * @throws DimensionMismatchException if derivative structures are inconsistent
87 * @throws NullArgumentException if x is null
88 */
89 @SafeVarargs
90 public final void addSamplePoint(final T x, final T[] ... value)
91 throws ZeroException, MathArithmeticException,
92 DimensionMismatchException, NullArgumentException {
93
94 NullArgumentException.check(x);
95 T factorial = x.getField().getOne();
96 for (int i = 0; i < value.length; ++i) {
97
98 final T[] y = value[i].clone();
99 if (i > 1) {
100 factorial = factorial.multiply(i);
101 final T inv = factorial.reciprocal();
102 for (int j = 0; j < y.length; ++j) {
103 y[j] = y[j].multiply(inv);
104 }
105 }
106
107 // update the bottom diagonal of the divided differences array
108 final int n = abscissae.size();
109 bottomDiagonal.add(n - i, y);
110 T[] bottom0 = y;
111 for (int j = i; j < n; ++j) {
112 final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
113 if (x.equals(abscissae.get(n - (j + 1)))) {
114 throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
115 }
116 final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal();
117 for (int k = 0; k < y.length; ++k) {
118 bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
119 }
120 bottom0 = bottom1;
121 }
122
123 // update the top diagonal of the divided differences array
124 topDiagonal.add(bottom0.clone());
125
126 // update the abscissae array
127 abscissae.add(x);
128 }
129 }
130
131 /** Interpolate value at a specified abscissa.
132 * @param x interpolation abscissa
133 * @return interpolated value
134 * @exception NoDataException if sample is empty
135 * @throws NullArgumentException if x is null
136 */
137 public T[] value(T x) throws NoDataException, NullArgumentException {
138
139 // safety check
140 NullArgumentException.check(x);
141 if (abscissae.isEmpty()) {
142 throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
143 }
144
145 final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
146 T valueCoeff = x.getField().getOne();
147 for (int i = 0; i < topDiagonal.size(); ++i) {
148 T[] dividedDifference = topDiagonal.get(i);
149 for (int k = 0; k < value.length; ++k) {
150 value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff));
151 }
152 final T deltaX = x.subtract(abscissae.get(i));
153 valueCoeff = valueCoeff.multiply(deltaX);
154 }
155
156 return value;
157 }
158
159 /** Interpolate value and first derivatives at a specified abscissa.
160 * @param x interpolation abscissa
161 * @param order maximum derivation order
162 * @return interpolated value and derivatives (value in row 0,
163 * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
164 * @exception NoDataException if sample is empty
165 * @throws NullArgumentException if x is null
166 */
167 public T[][] derivatives(T x, int order) throws NoDataException, NullArgumentException {
168
169 // safety check
170 NullArgumentException.check(x);
171 if (abscissae.isEmpty()) {
172 throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
173 }
174
175 final T zero = x.getField().getZero();
176 final T one = x.getField().getOne();
177 final T[] tj = MathArrays.buildArray(x.getField(), order + 1);
178 tj[0] = zero;
179 for (int i = 0; i < order; ++i) {
180 tj[i + 1] = tj[i].add(one);
181 }
182
183 final T[][] derivatives =
184 MathArrays.buildArray(x.getField(), order + 1, topDiagonal.get(0).length);
185 final T[] valueCoeff = MathArrays.buildArray(x.getField(), order + 1);
186 valueCoeff[0] = x.getField().getOne();
187 for (int i = 0; i < topDiagonal.size(); ++i) {
188 T[] dividedDifference = topDiagonal.get(i);
189 final T deltaX = x.subtract(abscissae.get(i));
190 for (int j = order; j >= 0; --j) {
191 for (int k = 0; k < derivatives[j].length; ++k) {
192 derivatives[j][k] =
193 derivatives[j][k].add(dividedDifference[k].multiply(valueCoeff[j]));
194 }
195 valueCoeff[j] = valueCoeff[j].multiply(deltaX);
196 if (j > 0) {
197 valueCoeff[j] = valueCoeff[j].add(tj[j].multiply(valueCoeff[j - 1]));
198 }
199 }
200 }
201
202 return derivatives;
203 }
204 }