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 }