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 }