View Javadoc
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.polynomials;
23  
24  import org.hipparchus.analysis.UnivariateFunction;
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.util.FastMath;
28  import org.hipparchus.util.MathArrays;
29  
30  /**
31   * Implements the representation of a real polynomial function in
32   * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
33   * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
34   * Analysis</b>, ISBN 038795452X, chapter 2.
35   * <p>
36   * The approximated function should be smooth enough for Lagrange polynomial
37   * to work well. Otherwise, consider using splines instead.</p>
38   *
39   */
40  public class PolynomialFunctionLagrangeForm implements UnivariateFunction {
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 double coefficients[];
47      /**
48       * Interpolating points (abscissas).
49       */
50      private final double x[];
51      /**
52       * Function values at interpolating points.
53       */
54      private final double 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 are not 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
71       * if two abscissae have the same value.
72       */
73      public PolynomialFunctionLagrangeForm(double x[], double y[])
74          throws MathIllegalArgumentException {
75          this.x = new double[x.length];
76          this.y = new double[y.length];
77          System.arraycopy(x, 0, this.x, 0, x.length);
78          System.arraycopy(y, 0, this.y, 0, y.length);
79          coefficientsComputed = false;
80  
81          if (!verifyInterpolationArray(x, y, false)) {
82              MathArrays.sortInPlace(this.x, this.y);
83              // Second check in case some abscissa is duplicated.
84              verifyInterpolationArray(this.x, this.y, true);
85          }
86      }
87  
88      /**
89       * Calculate the function value at the given point.
90       *
91       * @param z Point at which the function value is to be computed.
92       * @return the function value.
93       * @throws MathIllegalArgumentException if {@code x} and {@code y} have
94       * different lengths.
95       * @throws org.hipparchus.exception.MathIllegalArgumentException
96       * if {@code x} is not sorted in strictly increasing order.
97       * @throws MathIllegalArgumentException if the size of {@code x} is less
98       * than 2.
99       */
100     @Override
101     public double value(double z) {
102         return evaluateInternal(x, y, z);
103     }
104 
105     /**
106      * Returns the degree of the polynomial.
107      *
108      * @return the degree of the polynomial
109      */
110     public int degree() {
111         return x.length - 1;
112     }
113 
114     /**
115      * Returns a copy of the interpolating points array.
116      * <p>
117      * Changes made to the returned copy will not affect the polynomial.</p>
118      *
119      * @return a fresh copy of the interpolating points array
120      */
121     public double[] getInterpolatingPoints() {
122         double[] out = new double[x.length];
123         System.arraycopy(x, 0, out, 0, x.length);
124         return out;
125     }
126 
127     /**
128      * Returns a copy of the interpolating values array.
129      * <p>
130      * Changes made to the returned copy will not affect the polynomial.</p>
131      *
132      * @return a fresh copy of the interpolating values array
133      */
134     public double[] getInterpolatingValues() {
135         double[] out = new double[y.length];
136         System.arraycopy(y, 0, out, 0, y.length);
137         return out;
138     }
139 
140     /**
141      * Returns a copy of the coefficients array.
142      * <p>
143      * Changes made to the returned copy will not affect the polynomial.</p>
144      * <p>
145      * Note that coefficients computation can be ill-conditioned. Use with caution
146      * and only when it is necessary.</p>
147      *
148      * @return a fresh copy of the coefficients array
149      */
150     public double[] getCoefficients() {
151         if (!coefficientsComputed) {
152             computeCoefficients();
153         }
154         double[] out = new double[coefficients.length];
155         System.arraycopy(coefficients, 0, out, 0, coefficients.length);
156         return out;
157     }
158 
159     /**
160      * Evaluate the Lagrange polynomial using
161      * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
162      * Neville's Algorithm</a>. It takes O(n^2) time.
163      *
164      * @param x Interpolating points array.
165      * @param y Interpolating values array.
166      * @param z Point at which the function value is to be computed.
167      * @return the function value.
168      * @throws MathIllegalArgumentException if {@code x} and {@code y} have
169      * different lengths.
170      * @throws MathIllegalArgumentException
171      * if {@code x} is not sorted in strictly increasing order.
172      * @throws MathIllegalArgumentException if the size of {@code x} is less
173      * than 2.
174      */
175     public static double evaluate(double x[], double y[], double z)
176         throws MathIllegalArgumentException {
177         if (verifyInterpolationArray(x, y, false)) {
178             return evaluateInternal(x, y, z);
179         }
180 
181         // Array is not sorted.
182         final double[] xNew = new double[x.length];
183         final double[] yNew = new double[y.length];
184         System.arraycopy(x, 0, xNew, 0, x.length);
185         System.arraycopy(y, 0, yNew, 0, y.length);
186 
187         MathArrays.sortInPlace(xNew, yNew);
188         // Second check in case some abscissa is duplicated.
189         verifyInterpolationArray(xNew, yNew, true);
190         return evaluateInternal(xNew, yNew, z);
191     }
192 
193     /**
194      * Evaluate the Lagrange polynomial using
195      * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
196      * Neville's Algorithm</a>. It takes O(n^2) time.
197      *
198      * @param x Interpolating points array.
199      * @param y Interpolating values array.
200      * @param z Point at which the function value is to be computed.
201      * @return the function value.
202      * @throws MathIllegalArgumentException if {@code x} and {@code y} have
203      * different lengths.
204      * @throws org.hipparchus.exception.MathIllegalArgumentException
205      * if {@code x} is not sorted in strictly increasing order.
206      * @throws MathIllegalArgumentException if the size of {@code x} is less
207      * than 2.
208      */
209     private static double evaluateInternal(double x[], double y[], double z) {
210         int nearest = 0;
211         final int n = x.length;
212         final double[] c = new double[n];
213         final double[] d = new double[n];
214         double min_dist = Double.POSITIVE_INFINITY;
215         for (int i = 0; i < n; i++) {
216             // initialize the difference arrays
217             c[i] = y[i];
218             d[i] = y[i];
219             // find out the abscissa closest to z
220             final double dist = FastMath.abs(z - x[i]);
221             if (dist < min_dist) {
222                 nearest = i;
223                 min_dist = dist;
224             }
225         }
226 
227         // initial approximation to the function value at z
228         double value = y[nearest];
229 
230         for (int i = 1; i < n; i++) {
231             for (int j = 0; j < n-i; j++) {
232                 final double tc = x[j] - z;
233                 final double td = x[i+j] - z;
234                 final double divider = x[j] - x[i+j];
235                 // update the difference arrays
236                 final double w = (c[j+1] - d[j]) / divider;
237                 c[j] = tc * w;
238                 d[j] = td * w;
239             }
240             // sum up the difference terms to get the final value
241             if (nearest < 0.5*(n-i+1)) {
242                 value += c[nearest];    // fork down
243             } else {
244                 nearest--;
245                 value += d[nearest];    // fork up
246             }
247         }
248 
249         return value;
250     }
251 
252     /**
253      * Calculate the coefficients of Lagrange polynomial from the
254      * interpolation data. It takes O(n^2) time.
255      * Note that this computation can be ill-conditioned: Use with caution
256      * and only when it is necessary.
257      */
258     protected void computeCoefficients() {
259         final int n = degree() + 1;
260         coefficients = new double[n];
261         for (int i = 0; i < n; i++) {
262             coefficients[i] = 0.0;
263         }
264 
265         // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
266         final double[] c = new double[n+1];
267         c[0] = 1.0;
268         for (int i = 0; i < n; i++) {
269             for (int j = i; j > 0; j--) {
270                 c[j] = c[j-1] - c[j] * x[i];
271             }
272             c[0] *= -x[i];
273             c[i+1] = 1;
274         }
275 
276         final double[] tc = new double[n];
277         for (int i = 0; i < n; i++) {
278             // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
279             double d = 1;
280             for (int j = 0; j < n; j++) {
281                 if (i != j) {
282                     d *= x[i] - x[j];
283                 }
284             }
285             final double t = y[i] / d;
286             // Lagrange polynomial is the sum of n terms, each of which is a
287             // polynomial of degree n-1. tc[] are the coefficients of the i-th
288             // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
289             tc[n-1] = c[n];     // actually c[n] = 1
290             coefficients[n-1] += t * tc[n-1];
291             for (int j = n-2; j >= 0; j--) {
292                 tc[j] = c[j+1] + tc[j+1] * x[i];
293                 coefficients[j] += t * tc[j];
294             }
295         }
296 
297         coefficientsComputed = true;
298     }
299 
300     /**
301      * Check that the interpolation arrays are valid.
302      * The arrays features checked by this method are that both arrays have the
303      * same length and this length is at least 2.
304      *
305      * @param x Interpolating points array.
306      * @param y Interpolating values array.
307      * @param abort Whether to throw an exception if {@code x} is not sorted.
308      * @throws MathIllegalArgumentException if the array lengths are different.
309      * @throws MathIllegalArgumentException if the number of points is less than 2.
310      * @throws org.hipparchus.exception.MathIllegalArgumentException
311      * if {@code x} is not sorted in strictly increasing order and {@code abort}
312      * is {@code true}.
313      * @return {@code false} if the {@code x} is not sorted in increasing order,
314      * {@code true} otherwise.
315      * @see #evaluate(double[], double[], double)
316      * @see #computeCoefficients()
317      */
318     public static boolean verifyInterpolationArray(double x[], double y[], boolean abort)
319         throws MathIllegalArgumentException {
320         MathArrays.checkEqualLength(x, y);
321         if (x.length < 2) {
322             throw new MathIllegalArgumentException(LocalizedCoreFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
323         }
324 
325         return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
326     }
327 }