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    }