1 /*
2 * Licensed to the Hipparchus project 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 Hipparchus project 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 package org.hipparchus.analysis.polynomials;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
22 import org.hipparchus.exception.LocalizedCoreFormats;
23 import org.hipparchus.exception.MathIllegalArgumentException;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.MathArrays;
26
27 /**
28 * Implements the representation of a real polynomial function in
29 * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
30 * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
31 * Analysis</b>, ISBN 038795452X, chapter 2.
32 * <p>
33 * The approximated function should be smooth enough for Lagrange polynomial
34 * to work well. Otherwise, consider using splines instead.</p>
35 * @see PolynomialFunctionLagrangeForm
36 * @since 4.0
37 * @param <T> type of the field elements
38 */
39 public class FieldPolynomialFunctionLagrangeForm<T extends CalculusFieldElement<T>>
40 implements CalculusFieldUnivariateFunction<T> {
41 /**
42 * The coefficients of the polynomial, ordered by degree -- i.e.
43 * coefficients[0] is the constant term and coefficients[n] is the
44 * coefficient of x^n where n is the degree of the polynomial.
45 */
46 private T[] coefficients;
47 /**
48 * Interpolating points (abscissas).
49 */
50 private final T[] x;
51 /**
52 * Function values at interpolating points.
53 */
54 private final T[] y;
55 /**
56 * Whether the polynomial coefficients are available.
57 */
58 private boolean coefficientsComputed;
59
60 /**
61 * Construct a Lagrange polynomial with the given abscissas and function
62 * values. The order of interpolating points is important.
63 * <p>
64 * The constructor makes copy of the input arrays and assigns them.</p>
65 *
66 * @param x interpolating points
67 * @param y function values at interpolating points
68 * @throws MathIllegalArgumentException if the array lengths are different.
69 * @throws MathIllegalArgumentException if the number of points is less than 2.
70 * @throws MathIllegalArgumentException if two abscissae have the same value.
71 * @throws MathIllegalArgumentException if the abscissae are not sorted.
72 */
73 public FieldPolynomialFunctionLagrangeForm(final T[] x, final T[] y)
74 throws MathIllegalArgumentException {
75 this.x = x.clone();
76 this.y = y.clone();
77 coefficientsComputed = false;
78
79 MathArrays.checkEqualLength(x, y);
80 if (x.length < 2) {
81 throw new MathIllegalArgumentException(LocalizedCoreFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
82 }
83 MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, true);
84 }
85
86 /**
87 * Calculate the function value at the given point.
88 *
89 * @param z Point at which the function value is to be computed.
90 * @return the function value.
91 * @throws MathIllegalArgumentException if {@code x} and {@code y} have
92 * different lengths.
93 * @throws MathIllegalArgumentException
94 * if {@code x} is not sorted in strictly increasing order.
95 * @throws MathIllegalArgumentException if the size of {@code x} is less
96 * than 2.
97 */
98 @Override
99 public T value(final T z) {
100 int nearest = 0;
101 final int n = x.length;
102 final T[] c = y.clone();
103 final T[] d = c.clone();
104 double minDist = Double.POSITIVE_INFINITY;
105 for (int i = 0; i < n; i++) {
106 // find out the abscissa closest to z
107 final double dist = FastMath.abs(z.subtract(x[i])).getReal();
108 if (dist < minDist) {
109 nearest = i;
110 minDist = dist;
111 }
112 }
113
114 // initial approximation to the function value at z
115 T value = y[nearest];
116
117 for (int i = 1; i < n; i++) {
118 for (int j = 0; j < n-i; j++) {
119 final T tc = x[j].subtract(z);
120 final T td = x[i+j].subtract(z);
121 final T divider = x[j].subtract(x[i+j]);
122 // update the difference arrays
123 final T w = (c[j+1].subtract(d[j])).divide(divider);
124 c[j] = tc.multiply(w);
125 d[j] = td.multiply(w);
126 }
127 // sum up the difference terms to get the final value
128 if (nearest < 0.5*(n-i+1)) {
129 value = value.add(c[nearest]); // fork down
130 } else {
131 nearest--;
132 value = value.add(d[nearest]); // fork up
133 }
134 }
135
136 return value;
137 }
138
139 /**
140 * Returns the degree of the polynomial.
141 *
142 * @return the degree of the polynomial
143 */
144 public int degree() {
145 return x.length - 1;
146 }
147
148 /**
149 * Returns a copy of the interpolating points array.
150 * <p>
151 * Changes made to the returned copy will not affect the polynomial.</p>
152 *
153 * @return a fresh copy of the interpolating points array
154 */
155 public T[] getInterpolatingPoints() {
156 return x.clone();
157 }
158
159 /**
160 * Returns a copy of the interpolating values array.
161 * <p>
162 * Changes made to the returned copy will not affect the polynomial.</p>
163 *
164 * @return a fresh copy of the interpolating values array
165 */
166 public T[] getInterpolatingValues() {
167 return y.clone();
168 }
169
170 /**
171 * Returns a copy of the coefficients array.
172 * <p>
173 * Changes made to the returned copy will not affect the polynomial.</p>
174 * <p>
175 * Note that coefficients computation can be ill-conditioned. Use with caution
176 * and only when it is necessary.</p>
177 *
178 * @return a fresh copy of the coefficients array
179 */
180 public T[] getCoefficients() {
181 if (!coefficientsComputed) {
182 computeCoefficients();
183 }
184 return coefficients.clone();
185 }
186
187 /**
188 * Calculate the coefficients of Lagrange polynomial from the
189 * interpolation data. It takes O(n^2) time.
190 * Note that this computation can be ill-conditioned: Use with caution
191 * and only when it is necessary.
192 */
193 protected void computeCoefficients() {
194 final int n = degree() + 1;
195 final Field<T> field = x[0].getField();
196 coefficients = MathArrays.buildArray(field, n);
197
198 // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
199 final T[] c = MathArrays.buildArray(field, n + 1);
200 c[0] = field.getOne();
201 for (int i = 0; i < n; i++) {
202 for (int j = i; j > 0; j--) {
203 c[j] = c[j-1].subtract(c[j].multiply(x[i]));
204 }
205 c[0] = c[0].multiply(x[i].negate());
206 c[i+1] = field.getOne();
207 }
208
209 final T[] tc = MathArrays.buildArray(field, n);
210 for (int i = 0; i < n; i++) {
211 // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
212 T d = field.getOne();
213 for (int j = 0; j < n; j++) {
214 if (i != j) {
215 d = d.multiply(x[i].subtract(x[j]));
216 }
217 }
218 final T t = y[i].divide(d);
219 // Lagrange polynomial is the sum of n terms, each of which is a
220 // polynomial of degree n-1. tc[] are the coefficients of the i-th
221 // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
222 tc[n-1] = c[n]; // actually c[n] = 1
223 coefficients[n-1] = coefficients[n-1].add(t.multiply(tc[n-1]));
224 for (int j = n-2; j >= 0; j--) {
225 tc[j] = c[j+1].add(tc[j+1].multiply(x[i]));
226 coefficients[j] = coefficients[j].add(t.multiply(tc[j]));
227 }
228 }
229
230 coefficientsComputed = true;
231 }
232 }