View Javadoc

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    *      http://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.apache.commons.math.analysis.polynomials;
18  
19  import org.apache.commons.math.DuplicateSampleAbscissaException;
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MathRuntimeException;
22  import org.apache.commons.math.analysis.UnivariateRealFunction;
23  
24  /**
25   * Implements the representation of a real polynomial function in
26   * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
27   * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
28   * Analysis</b>, ISBN 038795452X, chapter 2.
29   * <p>
30   * The approximated function should be smooth enough for Lagrange polynomial
31   * to work well. Otherwise, consider using splines instead.</p>
32   *
33   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
34   * @since 1.2
35   */
36  public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
37  
38      /**
39       * The coefficients of the polynomial, ordered by degree -- i.e.
40       * coefficients[0] is the constant term and coefficients[n] is the 
41       * coefficient of x^n where n is the degree of the polynomial.
42       */
43      private double coefficients[];
44  
45      /**
46       * Interpolating points (abscissas) and the function values at these points.
47       */
48      private double x[], y[];
49  
50      /**
51       * Whether the polynomial coefficients are available.
52       */
53      private boolean coefficientsComputed;
54  
55      /**
56       * Construct a Lagrange polynomial with the given abscissas and function
57       * values. The order of interpolating points are not important.
58       * <p>
59       * The constructor makes copy of the input arrays and assigns them.</p>
60       * 
61       * @param x interpolating points
62       * @param y function values at interpolating points
63       * @throws IllegalArgumentException if input arrays are not valid
64       */
65      public PolynomialFunctionLagrangeForm(double x[], double y[])
66          throws IllegalArgumentException {
67  
68          verifyInterpolationArray(x, y);
69          this.x = new double[x.length];
70          this.y = new double[y.length];
71          System.arraycopy(x, 0, this.x, 0, x.length);
72          System.arraycopy(y, 0, this.y, 0, y.length);
73          coefficientsComputed = false;
74      }
75  
76      /**
77       * Calculate the function value at the given point.
78       *
79       * @param z the point at which the function value is to be computed
80       * @return the function value
81       * @throws FunctionEvaluationException if a runtime error occurs
82       * @see UnivariateRealFunction#value(double)
83       */
84      public double value(double z) throws FunctionEvaluationException {
85          try {
86              return evaluate(x, y, z);
87          } catch (DuplicateSampleAbscissaException e) {
88              throw new FunctionEvaluationException(e, z, e.getPattern(), e.getArguments());
89          }
90      }
91  
92      /**
93       * Returns the degree of the polynomial.
94       * 
95       * @return the degree of the polynomial
96       */
97      public int degree() {
98          return x.length - 1;
99      }
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      * 
132      * @return a fresh copy of the coefficients array
133      */
134     public double[] getCoefficients() {
135         if (!coefficientsComputed) {
136             computeCoefficients();
137         }
138         double[] out = new double[coefficients.length];
139         System.arraycopy(coefficients, 0, out, 0, coefficients.length);
140         return out;
141     }
142 
143     /**
144      * Evaluate the Lagrange polynomial using 
145      * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
146      * Neville's Algorithm</a>. It takes O(N^2) time.
147      * <p>
148      * This function is made public static so that users can call it directly
149      * without instantiating PolynomialFunctionLagrangeForm object.</p>
150      *
151      * @param x the interpolating points array
152      * @param y the interpolating values array
153      * @param z the point at which the function value is to be computed
154      * @return the function value
155      * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas
156      * @throws IllegalArgumentException if inputs are not valid
157      */
158     public static double evaluate(double x[], double y[], double z) throws
159         DuplicateSampleAbscissaException, IllegalArgumentException {
160 
161         int i, j, n, nearest = 0;
162         double value, c[], d[], tc, td, divider, w, dist, min_dist;
163 
164         verifyInterpolationArray(x, y);
165 
166         n = x.length;
167         c = new double[n];
168         d = new double[n];
169         min_dist = Double.POSITIVE_INFINITY;
170         for (i = 0; i < n; i++) {
171             // initialize the difference arrays
172             c[i] = y[i];
173             d[i] = y[i];
174             // find out the abscissa closest to z
175             dist = Math.abs(z - x[i]);
176             if (dist < min_dist) {
177                 nearest = i;
178                 min_dist = dist;
179             }
180         }
181 
182         // initial approximation to the function value at z
183         value = y[nearest];
184 
185         for (i = 1; i < n; i++) {
186             for (j = 0; j < n-i; j++) {
187                 tc = x[j] - z;
188                 td = x[i+j] - z;
189                 divider = x[j] - x[i+j];
190                 if (divider == 0.0) {
191                     // This happens only when two abscissas are identical.
192                     throw new DuplicateSampleAbscissaException(x[i], i, i+j);
193                 }
194                 // update the difference arrays
195                 w = (c[j+1] - d[j]) / divider;
196                 c[j] = tc * w;
197                 d[j] = td * w;
198             }
199             // sum up the difference terms to get the final value
200             if (nearest < 0.5*(n-i+1)) {
201                 value += c[nearest];    // fork down
202             } else {
203                 nearest--;
204                 value += d[nearest];    // fork up
205             }
206         }
207 
208         return value;
209     }
210 
211     /**
212      * Calculate the coefficients of Lagrange polynomial from the
213      * interpolation data. It takes O(N^2) time.
214      * <p>
215      * Note this computation can be ill-conditioned. Use with caution
216      * and only when it is necessary.</p>
217      *
218      * @throws ArithmeticException if any abscissas coincide
219      */
220     protected void computeCoefficients() throws ArithmeticException {
221         int i, j, n;
222         double c[], tc[], d, t;
223 
224         n = degree() + 1;
225         coefficients = new double[n];
226         for (i = 0; i < n; i++) {
227             coefficients[i] = 0.0;
228         }
229 
230         // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
231         c = new double[n+1];
232         c[0] = 1.0;
233         for (i = 0; i < n; i++) {
234             for (j = i; j > 0; j--) {
235                 c[j] = c[j-1] - c[j] * x[i];
236             }
237             c[0] *= (-x[i]);
238             c[i+1] = 1;
239         }
240 
241         tc = new double[n];
242         for (i = 0; i < n; i++) {
243             // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
244             d = 1;
245             for (j = 0; j < n; j++) {
246                 if (i != j) {
247                     d *= (x[i] - x[j]);
248                 }
249             }
250             if (d == 0.0) {
251                 // This happens only when two abscissas are identical.
252                 for (int k = 0; k < n; ++k) {
253                     if ((i != k) && (x[i] == x[k])) {
254                         throw MathRuntimeException.createArithmeticException("identical abscissas x[{0}] == x[{1}] == {2} cause division by zero",
255                                                                              i, k, x[i]);
256                     }
257                 }
258             }
259             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 (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 interpolating points must be distinct. However it is not
278      * verified here, it is checked in evaluate() and computeCoefficients().</p>
279      * 
280      * @param x the interpolating points array
281      * @param y the interpolating values array
282      * @throws IllegalArgumentException if not valid
283      * @see #evaluate(double[], double[], double)
284      * @see #computeCoefficients()
285      */
286     public static void verifyInterpolationArray(double x[], double y[]) throws
287         IllegalArgumentException {
288 
289         if (Math.min(x.length, y.length) < 2) {
290             throw MathRuntimeException.createIllegalArgumentException(
291                   "{0} points are required, got only {1}",
292                   2, Math.min(x.length, y.length));
293         }
294         if (x.length != y.length) {
295             throw MathRuntimeException.createIllegalArgumentException(
296                   "dimension mismatch {0} != {1}", x.length, y.length);
297         }
298     }
299 }