#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define myPI 3.1415927
/*================================================*/
/* this function returns Lagrange interpolation at
the point xx if x are the nodes and y are the values in the nodes
*/
float lagrange(float xx, float* x, float* y, int N)
{
float yy=0, tmp;
int k,j;
for(k=0;k<N;k++)
{
tmp=1;
for(j=0;j<N;j++) if(j!=k) tmp = tmp* (xx-x[j])/(x[k]-x[j]);
yy = yy + y[k]*tmp;
};
return yy;
};
/*================================================*/
/* the function to be interpolated */
float my_function(float x)
{ /* return this value f(x) = x^1.5 */
return x*sqrt(x);
}
/*================================================*/
int main(void)
{
int Npol=2; /* this is the power of the polynomial to use */
int N;
float *x; /* array of the chebyshev knots */
float *y; /* array of the f(x)-values in the vortecies */
float xx, yy;
float a=0, b=1.0; /* interval of interpolation */
float h=0.1; /* step */
int i;
float tmp,err,maxerr;
N=Npol+1; /* number of points is the power of the polynom + 1 */
x = malloc(sizeof(float)*N);
y = malloc(sizeof(float)*N);
/* define the nodes and the function values in them */
printf("function to interpolate:\nx-nodes y-values\n");
for(i=0;i<N;i++)
{ /*x[i]=a + (b-a)*(float)i/(float)(N-1);*/ /*equdistant nodes */
/* Chebyshev nodes */
x[i] =0.5*((b-a)*cos((float)(2*i+1)/(float)(2*N)*myPI)+a+b);
y[i]=my_function(x[i]); /* these are the function values */
printf("%f %f\n",x[i],y[i]);
};
printf("-----------------\n");
maxerr=0;
printf("Results:\n");
printf(" x \t L(x) \t f(x) \t err\n");
for(xx=a;xx<=b;xx=xx+h)
{
yy=lagrange(xx,x,y,N); /* Lagrange interpolation */
tmp = my_function(xx);
err=fabs(yy-tmp);
if(err>maxerr) maxerr=err;
printf("%f\t%f\t%f\t%f\n",xx,yy,tmp,err);
} ;
printf("Max error on [%f; %f] is %f\n",a,b,maxerr);
free(x); free(y);
return 0;
}
|
