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 * Estimates coefficients of a polynomial passing through provided set of x and y points.
20 * This implementation is more accurate than other implementations such as implementations of
21 * {@link com.irurueta.numerical.polynomials.estimators.PolynomialEstimator}.
22 * Those implementations are based on curve fitting a polynomial of a given degree based on a number
23 * of points that can be much larger than the number of points, which results in a polynomial that
24 * does not exactly pass through provided points, but provides a minimum square error result.
25 * On the other hand, this implementation builds a polynomial of order equal to the length of
26 * provided x and y points.
27 * Ths implementation is more accurate than {@link SimpleInterpolatingPolynomialEstimator}, at the
28 * expense of greater computational cost (N^3 of this method vs N^2 of the simple one).
29 */
30 public class AccurateInterpolatingPolynomialEstimator extends InterpolatingPolynomialEstimator {
31
32 /**
33 * Estimates polynomial coefficients from provided x and y points.
34 *
35 * @param xa x points the estimated polynomial passes through.
36 * @param ya y points the estimated polynomial passes through.
37 * @param cof instance where coefficients of estimated polynomial will be stored.
38 * @throws IllegalArgumentException if any of the provided values doesn't have the same length.
39 * @throws InterpolationException if interpolation fails for numerical reasons.
40 */
41 @Override
42 public void estimate(final double[] xa, final double[] ya, final double[] cof) throws InterpolationException {
43 int k;
44 int j;
45 int i;
46 double xmin;
47 final var n = xa.length;
48 final var x = new double[n];
49 final var y = new double[n];
50
51 for (j = 0; j < n; j++) {
52 x[j] = xa[j];
53 y[j] = ya[j];
54 }
55 for (j = 0; j < n; j++) {
56 final var interp = new PolynomialInterpolator(x, y, n - j, false);
57 // extrapolate to x = 0
58 cof[j] = interp.rawinterp(0, 0.);
59 xmin = 1.0e99;
60 k = -1;
61 for (i = 0; i < n - j; i++) {
62 // Find the remaining xi of smallest absolute value
63 if (Math.abs(x[i]) < xmin) {
64 xmin = Math.abs(x[i]);
65 k = i;
66 }
67 if (x[i] != 0.0) {
68 // (meanwhile reducing all the terms)
69 y[i] = (y[i] - cof[j]) / x[i];
70 }
71 }
72 // and eliminate it
73 for (i = k + 1; i < n - j; i++) {
74 y[i - 1] = y[i];
75 x[i - 1] = x[i];
76 }
77 }
78 }
79 }