Here you can find the source of interpolate(double xa[], double ya[], double x)
| Parameter | Description |
|---|---|
| xa | the array of known x-values |
| ya | the array of known y-values |
| x | the x value for which the y value will be computed |
public static double interpolate(double xa[], double ya[], double x)
//package com.java2s; /*/*w w w.j av a 2 s . c om*/ JOpenChart Java Charting Library and Toolkit Copyright (C) 2001 Sebastian M?ller http://jopenchart.sourceforge.net This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ChartUtilities.java Created on 21. September 2001, 17:42 */ public class Main { /** This function performs a polynomial interpolation using a set of * given x and y values. It uses Neville's interpolation algorithm. * @param xa the array of known x-values * @param ya the array of known y-values * @param x the x value for which the y value will be computed * @return the corresponding y value */ public static double interpolate(double xa[], double ya[], double x) { /* Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y. If P(x) is the polynomial of degree N ? 1 such that P(xa[i]) = ya[i]; i = 1...n, then the returned value y = P(x). */ if (xa.length != ya.length || xa.length == 0 || ya.length == 0) { System.out.println("** Invalid Parameter"); return Double.NaN; } int n = xa.length; double y = 0.0; double dy = 0.0; int i, m, ns = 1; double den, dif, dift, ho, hp, w; double[] c = new double[n]; double[] d = new double[n]; dif = Math.abs(x - xa[0]); for (i = 0; i < n; i++) { // Here we find the index ns of the closest table entry, if ((dift = Math.abs(x - xa[i])) < dif) { ns = i; dif = dift; } c[i] = ya[i]; // and initialize the tableau of c's and d's. d[i] = ya[i]; } y = ya[ns--]; // This is the initial approximation to y. //System.out.println("** y ~ "+y); for (m = 0; m < n - 1; m++) { // For each column of the tableau, for (i = 0; i < n - m - 1; i++) { // we loop over the current c's and d's and update them. //System.out.println("** m = "+m+", i = "+i); ho = xa[i] - x; hp = xa[i + m + 1] - x; w = c[i + 1] - d[i]; if ((den = ho - hp) == 0.0) { return Double.NaN; } // This error can occur only if two input xa's are (to within roundof identical. //System.out.println("** ho = "+ho+", hp = "+hp); den = w / den; d[i] = hp * den; // Here the c's and d's are updated. c[i] = ho * den; //System.out.println("** c[i] = "+c[i]+", d[i] = "+d[i]); } y += (dy = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--])); //System.out.println("** dy = "+dy+", y = "+y); /* After each column in the tableau is completed, we decide which correction, c or d, we want to add to our accumulating value of y, i.e., which path to take through the tableau forking up or down. We do this in such a way as to take the most "straight line" route through the tableau to its apex, updating ns accordingly to keep track of where we are. This route keeps the partial approximations centered (insofar as possible) on the target x. The last dy added is thus the error indication. */ } return y; } }