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.math3.analysis.polynomials;
018
019 import org.apache.commons.math3.analysis.UnivariateFunction;
020 import org.apache.commons.math3.util.FastMath;
021 import org.apache.commons.math3.util.MathArrays;
022 import org.apache.commons.math3.exception.DimensionMismatchException;
023 import org.apache.commons.math3.exception.NumberIsTooSmallException;
024 import org.apache.commons.math3.exception.util.LocalizedFormats;
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 $Id: PolynomialFunctionLagrangeForm.java 1364387 2012-07-22 18:14:11Z tn $
036 * @since 1.2
037 */
038 public class PolynomialFunctionLagrangeForm implements UnivariateFunction {
039 /**
040 * The coefficients of the polynomial, ordered by degree -- i.e.
041 * coefficients[0] is the constant term and coefficients[n] is the
042 * coefficient of x^n where n is the degree of the polynomial.
043 */
044 private double coefficients[];
045 /**
046 * Interpolating points (abscissas).
047 */
048 private final double x[];
049 /**
050 * Function values at interpolating points.
051 */
052 private final double y[];
053 /**
054 * Whether the polynomial coefficients are available.
055 */
056 private boolean coefficientsComputed;
057
058 /**
059 * Construct a Lagrange polynomial with the given abscissas and function
060 * values. The order of interpolating points are not important.
061 * <p>
062 * The constructor makes copy of the input arrays and assigns them.</p>
063 *
064 * @param x interpolating points
065 * @param y function values at interpolating points
066 * @throws DimensionMismatchException if the array lengths are different.
067 * @throws NumberIsTooSmallException if the number of points is less than 2.
068 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
069 * if two abscissae have the same value.
070 */
071 public PolynomialFunctionLagrangeForm(double x[], double y[]) {
072 this.x = new double[x.length];
073 this.y = new double[y.length];
074 System.arraycopy(x, 0, this.x, 0, x.length);
075 System.arraycopy(y, 0, this.y, 0, y.length);
076 coefficientsComputed = false;
077
078 if (!verifyInterpolationArray(x, y, false)) {
079 MathArrays.sortInPlace(this.x, this.y);
080 // Second check in case some abscissa is duplicated.
081 verifyInterpolationArray(this.x, this.y, true);
082 }
083 }
084
085 /**
086 * Calculate the function value at the given point.
087 *
088 * @param z Point at which the function value is to be computed.
089 * @return the function value.
090 * @throws DimensionMismatchException if {@code x} and {@code y} have
091 * different lengths.
092 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
093 * if {@code x} is not sorted in strictly increasing order.
094 * @throws NumberIsTooSmallException if the size of {@code x} is less
095 * than 2.
096 */
097 public double value(double z) {
098 return evaluateInternal(x, y, z);
099 }
100
101 /**
102 * Returns the degree of the polynomial.
103 *
104 * @return the degree of the polynomial
105 */
106 public int degree() {
107 return x.length - 1;
108 }
109
110 /**
111 * Returns a copy of the interpolating points array.
112 * <p>
113 * Changes made to the returned copy will not affect the polynomial.</p>
114 *
115 * @return a fresh copy of the interpolating points array
116 */
117 public double[] getInterpolatingPoints() {
118 double[] out = new double[x.length];
119 System.arraycopy(x, 0, out, 0, x.length);
120 return out;
121 }
122
123 /**
124 * Returns a copy of the interpolating values array.
125 * <p>
126 * Changes made to the returned copy will not affect the polynomial.</p>
127 *
128 * @return a fresh copy of the interpolating values array
129 */
130 public double[] getInterpolatingValues() {
131 double[] out = new double[y.length];
132 System.arraycopy(y, 0, out, 0, y.length);
133 return out;
134 }
135
136 /**
137 * Returns a copy of the coefficients array.
138 * <p>
139 * Changes made to the returned copy will not affect the polynomial.</p>
140 * <p>
141 * Note that coefficients computation can be ill-conditioned. Use with caution
142 * and only when it is necessary.</p>
143 *
144 * @return a fresh copy of the coefficients array
145 */
146 public double[] getCoefficients() {
147 if (!coefficientsComputed) {
148 computeCoefficients();
149 }
150 double[] out = new double[coefficients.length];
151 System.arraycopy(coefficients, 0, out, 0, coefficients.length);
152 return out;
153 }
154
155 /**
156 * Evaluate the Lagrange polynomial using
157 * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
158 * Neville's Algorithm</a>. It takes O(n^2) time.
159 *
160 * @param x Interpolating points array.
161 * @param y Interpolating values array.
162 * @param z Point at which the function value is to be computed.
163 * @return the function value.
164 * @throws DimensionMismatchException if {@code x} and {@code y} have
165 * different lengths.
166 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
167 * if {@code x} is not sorted in strictly increasing order.
168 * @throws NumberIsTooSmallException if the size of {@code x} is less
169 * than 2.
170 */
171 public static double evaluate(double x[], double y[], double z) {
172 if (verifyInterpolationArray(x, y, false)) {
173 return evaluateInternal(x, y, z);
174 }
175
176 // Array is not sorted.
177 final double[] xNew = new double[x.length];
178 final double[] yNew = new double[y.length];
179 System.arraycopy(x, 0, xNew, 0, x.length);
180 System.arraycopy(y, 0, yNew, 0, y.length);
181
182 MathArrays.sortInPlace(xNew, yNew);
183 // Second check in case some abscissa is duplicated.
184 verifyInterpolationArray(xNew, yNew, true);
185 return evaluateInternal(xNew, yNew, z);
186 }
187
188 /**
189 * Evaluate the Lagrange polynomial using
190 * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
191 * Neville's Algorithm</a>. It takes O(n^2) time.
192 *
193 * @param x Interpolating points array.
194 * @param y Interpolating values array.
195 * @param z Point at which the function value is to be computed.
196 * @return the function value.
197 * @throws DimensionMismatchException if {@code x} and {@code y} have
198 * different lengths.
199 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
200 * if {@code x} is not sorted in strictly increasing order.
201 * @throws NumberIsTooSmallException if the size of {@code x} is less
202 * than 2.
203 */
204 private static double evaluateInternal(double x[], double y[], double z) {
205 int nearest = 0;
206 final int n = x.length;
207 final double[] c = new double[n];
208 final double[] d = new double[n];
209 double min_dist = Double.POSITIVE_INFINITY;
210 for (int i = 0; i < n; i++) {
211 // initialize the difference arrays
212 c[i] = y[i];
213 d[i] = y[i];
214 // find out the abscissa closest to z
215 final double dist = FastMath.abs(z - x[i]);
216 if (dist < min_dist) {
217 nearest = i;
218 min_dist = dist;
219 }
220 }
221
222 // initial approximation to the function value at z
223 double value = y[nearest];
224
225 for (int i = 1; i < n; i++) {
226 for (int j = 0; j < n-i; j++) {
227 final double tc = x[j] - z;
228 final double td = x[i+j] - z;
229 final double divider = x[j] - x[i+j];
230 // update the difference arrays
231 final double w = (c[j+1] - d[j]) / divider;
232 c[j] = tc * w;
233 d[j] = td * w;
234 }
235 // sum up the difference terms to get the final value
236 if (nearest < 0.5*(n-i+1)) {
237 value += c[nearest]; // fork down
238 } else {
239 nearest--;
240 value += d[nearest]; // fork up
241 }
242 }
243
244 return value;
245 }
246
247 /**
248 * Calculate the coefficients of Lagrange polynomial from the
249 * interpolation data. It takes O(n^2) time.
250 * Note that this computation can be ill-conditioned: Use with caution
251 * and only when it is necessary.
252 */
253 protected void computeCoefficients() {
254 final int n = degree() + 1;
255 coefficients = new double[n];
256 for (int i = 0; i < n; i++) {
257 coefficients[i] = 0.0;
258 }
259
260 // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
261 final double[] c = new double[n+1];
262 c[0] = 1.0;
263 for (int i = 0; i < n; i++) {
264 for (int j = i; j > 0; j--) {
265 c[j] = c[j-1] - c[j] * x[i];
266 }
267 c[0] *= -x[i];
268 c[i+1] = 1;
269 }
270
271 final double[] tc = new double[n];
272 for (int i = 0; i < n; i++) {
273 // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
274 double d = 1;
275 for (int j = 0; j < n; j++) {
276 if (i != j) {
277 d *= x[i] - x[j];
278 }
279 }
280 final double t = y[i] / d;
281 // Lagrange polynomial is the sum of n terms, each of which is a
282 // polynomial of degree n-1. tc[] are the coefficients of the i-th
283 // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
284 tc[n-1] = c[n]; // actually c[n] = 1
285 coefficients[n-1] += t * tc[n-1];
286 for (int j = n-2; j >= 0; j--) {
287 tc[j] = c[j+1] + tc[j+1] * x[i];
288 coefficients[j] += t * tc[j];
289 }
290 }
291
292 coefficientsComputed = true;
293 }
294
295 /**
296 * Check that the interpolation arrays are valid.
297 * The arrays features checked by this method are that both arrays have the
298 * same length and this length is at least 2.
299 *
300 * @param x Interpolating points array.
301 * @param y Interpolating values array.
302 * @param abort Whether to throw an exception if {@code x} is not sorted.
303 * @throws DimensionMismatchException if the array lengths are different.
304 * @throws NumberIsTooSmallException if the number of points is less than 2.
305 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
306 * if {@code x} is not sorted in strictly increasing order and {@code abort}
307 * is {@code true}.
308 * @return {@code false} if the {@code x} is not sorted in increasing order,
309 * {@code true} otherwise.
310 * @see #evaluate(double[], double[], double)
311 * @see #computeCoefficients()
312 */
313 public static boolean verifyInterpolationArray(double x[], double y[], boolean abort) {
314 if (x.length != y.length) {
315 throw new DimensionMismatchException(x.length, y.length);
316 }
317 if (x.length < 2) {
318 throw new NumberIsTooSmallException(LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
319 }
320
321 return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
322 }
323 }