    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.analysis.polynomials;
    018
    019    import org.apache.commons.math.DuplicateSampleAbscissaException;
    020    import org.apache.commons.math.MathRuntimeException;
    021    import org.apache.commons.math.analysis.UnivariateRealFunction;
    022    import org.apache.commons.math.FunctionEvaluationException;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.util.FastMath;
    025
    026    /**
    027     * Implements the representation of a real polynomial function in
    028     * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
    029     * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
    030     * Analysis</b>, ISBN 038795452X, chapter 2.
    031     * <p>
    032     * The approximated function should be smooth enough for Lagrange polynomial
    033     * to work well. Otherwise, consider using splines instead.</p>
    034     *
    035     * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 f??vr. 2011) $
    036     * @since 1.2
    037     */
    038    public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
    039
    040        /**
    041         * The coefficients of the polynomial, ordered by degree -- i.e.
    042         * coefficients[0] is the constant term and coefficients[n] is the
    043         * coefficient of x^n where n is the degree of the polynomial.
    044         */
    045        private double coefficients[];
    046
    047        /**
    048         * Interpolating points (abscissas).
    049         */
    050        private final double x[];
    051
    052        /**
    053         * Function values at interpolating points.
    054         */
    055        private final double y[];
    056
    057        /**
    058         * Whether the polynomial coefficients are available.
    059         */
    060        private boolean coefficientsComputed;
    061
    062        /**
    063         * Construct a Lagrange polynomial with the given abscissas and function
    064         * values. The order of interpolating points are not important.
    065         * <p>
    066         * The constructor makes copy of the input arrays and assigns them.</p>
    067         *
    068         * @param x interpolating points
    069         * @param y function values at interpolating points
    070         * @throws IllegalArgumentException if input arrays are not valid
    071         */
    072        public PolynomialFunctionLagrangeForm(double x[], double y[])
    073            throws IllegalArgumentException {
    074
    075            verifyInterpolationArray(x, y);
    076            this.x = new double[x.length];
    077            this.y = new double[y.length];
    078            System.arraycopy(x, 0, this.x, 0, x.length);
    079            System.arraycopy(y, 0, this.y, 0, y.length);
    080            coefficientsComputed = false;
    081        }
    082
    083        /** {@inheritDoc} */
    084        public double value(double z) throws FunctionEvaluationException {
    085            try {
    086                return evaluate(x, y, z);
    087            } catch (DuplicateSampleAbscissaException e) {
    088                throw new FunctionEvaluationException(z, e.getSpecificPattern(), e.getGeneralPattern(), e.getArguments());
    089            }
    090        }
    091
    092        /**
    093         * Returns the degree of the polynomial.
    094         *
    095         * @return the degree of the polynomial
    096         */
    097        public int degree() {
    098            return x.length - 1;
    099        }
    100
    101        /**
    102         * Returns a copy of the interpolating points array.
    103         * <p>
    104         * Changes made to the returned copy will not affect the polynomial.</p>
    105         *
    106         * @return a fresh copy of the interpolating points array
    107         */
    108        public double[] getInterpolatingPoints() {
    109            double[] out = new double[x.length];
    110            System.arraycopy(x, 0, out, 0, x.length);
    111            return out;
    112        }
    113
    114        /**
    115         * Returns a copy of the interpolating values 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 values array
    120         */
    121        public double[] getInterpolatingValues() {
    122            double[] out = new double[y.length];
    123            System.arraycopy(y, 0, out, 0, y.length);
    124            return out;
    125        }
    126
    127        /**
    128         * Returns a copy of the coefficients array.
    129         * <p>
    130         * Changes made to the returned copy will not affect the polynomial.</p>
    131         * <p>
    132         * Note that coefficients computation can be ill-conditioned. Use with caution
    133         * and only when it is necessary.</p>
    134         *
    135         * @return a fresh copy of the coefficients array
    136         */
    137        public double[] getCoefficients() {
    138            if (!coefficientsComputed) {
    139                computeCoefficients();
    140            }
    141            double[] out = new double[coefficients.length];
    142            System.arraycopy(coefficients, 0, out, 0, coefficients.length);
    143            return out;
    144        }
    145
    146        /**
    147         * Evaluate the Lagrange polynomial using
    148         * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
    149         * Neville's Algorithm</a>. It takes O(N^2) time.
    150         * <p>
    151         * This function is made public static so that users can call it directly
    152         * without instantiating PolynomialFunctionLagrangeForm object.</p>
    153         *
    154         * @param x the interpolating points array
    155         * @param y the interpolating values array
    156         * @param z the point at which the function value is to be computed
    157         * @return the function value
    158         * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas
    159         * @throws IllegalArgumentException if inputs are not valid
    160         */
    161        public static double evaluate(double x[], double y[], double z) throws
    162            DuplicateSampleAbscissaException, IllegalArgumentException {
    163
    164            verifyInterpolationArray(x, y);
    165
    166            int nearest = 0;
    167            final int n = x.length;
    168            final double[] c = new double[n];
    169            final double[] d = new double[n];
    170            double min_dist = Double.POSITIVE_INFINITY;
    171            for (int i = 0; i < n; i++) {
    172                // initialize the difference arrays
    173                c[i] = y[i];
    174                d[i] = y[i];
    175                // find out the abscissa closest to z
    176                final double dist = FastMath.abs(z - x[i]);
    177                if (dist < min_dist) {
    178                    nearest = i;
    179                    min_dist = dist;
    180                }
    181            }
    182
    183            // initial approximation to the function value at z
    184            double value = y[nearest];
    185
    186            for (int i = 1; i < n; i++) {
    187                for (int j = 0; j < n-i; j++) {
    188                    final double tc = x[j] - z;
    189                    final double td = x[i+j] - z;
    190                    final double divider = x[j] - x[i+j];
    191                    if (divider == 0.0) {
    192                        // This happens only when two abscissas are identical.
    193                        throw new DuplicateSampleAbscissaException(x[i], i, i+j);
    194                    }
    195                    // update the difference arrays
    196                    final double w = (c[j+1] - d[j]) / divider;
    197                    c[j] = tc * w;
    198                    d[j] = td * w;
    199                }
    200                // sum up the difference terms to get the final value
    201                if (nearest < 0.5*(n-i+1)) {
    202                    value += c[nearest];    // fork down
    203                } else {
    204                    nearest--;
    205                    value += d[nearest];    // fork up
    206                }
    207            }
    208
    209            return value;
    210        }
    211
    212        /**
    213         * Calculate the coefficients of Lagrange polynomial from the
    214         * interpolation data. It takes O(N^2) time.
    215         * <p>
    216         * Note this computation can be ill-conditioned. Use with caution
    217         * and only when it is necessary.</p>
    218         *
    219         * @throws ArithmeticException if any abscissas coincide
    220         */
    221        protected void computeCoefficients() throws ArithmeticException {
    222
    223            final int n = degree() + 1;
    224            coefficients = new double[n];
    225            for (int i = 0; i < n; i++) {
    226                coefficients[i] = 0.0;
    227            }
    228
    229            // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
    230            final double[] c = new double[n+1];
    231            c[0] = 1.0;
    232            for (int i = 0; i < n; i++) {
    233                for (int j = i; j > 0; j--) {
    234                    c[j] = c[j-1] - c[j] * x[i];
    235                }
    236                c[0] *= -x[i];
    237                c[i+1] = 1;
    238            }
    239
    240            final double[] tc = new double[n];
    241            for (int i = 0; i < n; i++) {
    242                // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
    243                double d = 1;
    244                for (int j = 0; j < n; j++) {
    245                    if (i != j) {
    246                        d *= x[i] - x[j];
    247                    }
    248                }
    249                if (d == 0.0) {
    250                    // This happens only when two abscissas are identical.
    251                    for (int k = 0; k < n; ++k) {
    252                        if ((i != k) && (x[i] == x[k])) {
    253                            throw MathRuntimeException.createArithmeticException(
    254                                  LocalizedFormats.IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO,
    255                                  i, k, x[i]);
    256                        }
    257                    }
    258                }
    259                final double t = y[i] / d;
    260                // Lagrange polynomial is the sum of n terms, each of which is a
    261                // polynomial of degree n-1. tc[] are the coefficients of the i-th
    262                // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
    263                tc[n-1] = c[n];     // actually c[n] = 1
    264                coefficients[n-1] += t * tc[n-1];
    265                for (int j = n-2; j >= 0; j--) {
    266                    tc[j] = c[j+1] + tc[j+1] * x[i];
    267                    coefficients[j] += t * tc[j];
    268                }
    269            }
    270
    271            coefficientsComputed = true;
    272        }
    273
    274        /**
    275         * Verifies that the interpolation arrays are valid.
    276         * <p>
    277         * The arrays features checked by this method are that both arrays have the
    278         * same length and this length is at least 2.
    279         * </p>
    280         * <p>
    281         * The interpolating points must be distinct. However it is not
    282         * verified here, it is checked in evaluate() and computeCoefficients().
    283         * </p>
    284         *
    285         * @param x the interpolating points array
    286         * @param y the interpolating values array
    287         * @throws IllegalArgumentException if not valid
    288         * @see #evaluate(double[], double[], double)
    289         * @see #computeCoefficients()
    290         */
    291        public static void verifyInterpolationArray(double x[], double y[])
    292            throws IllegalArgumentException {
    293
    294            if (x.length != y.length) {
    295                throw MathRuntimeException.createIllegalArgumentException(
    296                      LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, x.length, y.length);
    297            }
    298
    299            if (x.length < 2) {
    300                throw MathRuntimeException.createIllegalArgumentException(
    301                      LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length);
    302            }
    303
    304        }
    305    }