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