interpolation.cpp

Neville插值算法

#include <stdio.h>
#include <math.h>
#include <malloc.h>

double polint(double *xa, double *ya, int n, double x);

double polint(double *xa, double *ya, int n, double x)
{
	int i, m, ns;
	double den, dif, dift, ho, hp, w, y;
	double *c, *d;

	ns = 1;
	dif = fabs(x-xa[1]);

	c = (double*)malloc(sizeof(double)*(n+1));
	d = (double*)malloc(sizeof(double)*(n+1));

	for(i=1; i<=n; i++){
		if((dift=fabs(x-xa[i]))<dif){
			ns = i; 
			dif = dift;
		}
		c[i] = ya[i];
		d[i] = ya[i];
	}

	y = ya[ns--];
	
	for(m=1; m<n; m++){
		for(i=1; i<=n-m; i++){
			ho = xa[i]-x;
			hp = xa[i+m]-x;
			w = c[i+1]-d[i];
			if((den=ho-hp)==0.0){
				printf("Error in routine polint\n");
			}
			den = w/den;
			d[i] = hp*den;
			c[i] = ho*den;
		}
		y += (2*ns<(n-m)?c[ns+1]:d[ns--]);
	}
	free(d);
	free(c);

	return y;
}

void main()
{
	int i, n;
	double x, y, y_polin;
	double xa[5] = {0};
	double ya[5] = {0};

	n = 5;

	for(i=1; i<=5; i++){
		xa[i] = i;
		ya[i] = xa[i]*xa[i] - 3.5*xa[i] + 10;

	}

	x = 3.95;
	y = x*x - 3.5*x + 10;
	y_polin = polint(xa, ya, n, x);
	printf("y = %f\n", y);
	printf("y_polin = %f\n", y_polin);
}