This function performs a polynomial interpolation using a set of given x and y values. - Java java.lang

Java examples for java.lang:Math Algorithm

Description

This function performs a polynomial interpolation using a set of given x and y values.

Demo Code

/*/*from w w  w. j  a va2s. c  o  m*/
    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
 */
import java.util.*;

public class Main{
    public static void main(String[] argv) throws Exception{
        double xa = 2.45678;
        double ya = 2.45678;
        double x = 2.45678;
        System.out.println(interpolate(xa,ya,x));
    }
    /** 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;
    }
}

Related Tutorials