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