FieldPolynomialFunctionLagrangeForm.java

  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;

  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;

  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;

  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;

  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.     }

  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.         }

  114.         // initial approximation to the function value at z
  115.         T value = y[nearest];

  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.         }

  136.         return value;
  137.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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);

  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.         }

  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.         }

  230.         coefficientsComputed = true;
  231.     }
  232. }
Created with JaCoCo 0.8.12.202403310830