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 }