PolynomialFunctionLagrangeForm.java

  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.  */

  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;

  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;

  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;

  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;

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

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

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

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

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

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

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

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

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

  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 minDist = 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 < minDist) {
  222.                 nearest = i;
  223.                 minDist = dist;
  224.             }
  225.         }

  227.         // initial approximation to the function value at z
  228.         double value = y[nearest];

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

  249.         return value;
  250.     }

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

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

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

  297.         coefficientsComputed = true;
  298.     }

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

  325.         return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
  326.     }
  327. }
Created with JaCoCo 0.8.12.202403310830