Newton's Forward Interpolation : Numerical Analysis C Program
#include<stdio.h>
#include<malloc.h>
#include<stdlib.h>
void differenceTable(double **A,double **L,int ROW);
double result(double **A,double **L,double x,int ROW);
double factorial(int n);
FILE *fp;
int main(int argc, char *argv[]) {
double **A,*b,*xn,*xnplus1,*p,temp,x,**L,**U,lambda,error,t1=0.0;
int ROW,i,j,k,t=1,m=1,bo=1;
char ch;
if(argc==1) {
printf("\n Enter the number of values:");
scanf("%d",&ROW);
printf("\n Enter the value at which value has to be aprroximated:");
scanf("%lf",&x);
A=(double **)malloc((ROW+1)*sizeof(double*));
L=(double **)malloc((ROW+1)*sizeof(double*));
U=(double **)malloc((ROW+1)*sizeof(double*));
for(i=0;i<=ROW;i++) {
A[i]=(double *)malloc((ROW+1)*sizeof(double));
L[i]=(double *)malloc((ROW+1)*sizeof(double));
U[i]=(double *)malloc((ROW+1)*sizeof(double));
}
b=(double *)malloc((ROW+1)*sizeof(double));
p=(double *)malloc((ROW+1)*sizeof(double));
xn=(double *)malloc((ROW+1)*sizeof(double));
xnplus1=(double *)malloc((ROW+1)*sizeof(double));
for(i=1;i<=ROW;i++) {
printf("\n Enter the value of X[%d]:",i);
scanf("%lf",&A[i][1]);
}
for(i=1;i<=ROW;i++) {
printf("\n Enter the value of Y[%d]:",i);
scanf("%lf",&A[i][2]);
}
}
else if (argc==2) {
fp=fopen(argv[1],"r");
if(fp==NULL) {
printf("\n File not found, program will terminate:");
exit(0);
}
fscanf(fp,"%d",&ROW);
fscanf(fp,"%lf",&x);
A=(double **)malloc((ROW+1)*sizeof(double*));
L=(double **)malloc((ROW+1)*sizeof(double*));
U=(double **)malloc((ROW+1)*sizeof(double*));
for(i=0;i<=ROW;i++) {
A[i]=(double *)malloc((ROW+1)*sizeof(double));
L[i]=(double *)malloc((ROW+1)*sizeof(double));
U[i]=(double *)malloc((ROW+1)*sizeof(double));
}
b=(double *)malloc((ROW+1)*sizeof(double));
p=(double *)malloc((ROW+1)*sizeof(double));
xn=(double *)malloc((ROW+1)*sizeof(double));
xnplus1=(double *)malloc((ROW+1)*sizeof(double));
while(!feof(fp)) {
for(i=1;i<=ROW;i++) {
for(j=1;j<=2;j++) {
fscanf(fp,"%lf ",&A[i][j]);
}
}
}
fclose(fp);
}
else {
printf("\n Invalid Arguments,program will terminate:\n");
exit(0);
}
printf("\n X Y");
printf("\n------------------------------------------------------------------\n");
for(i=1;i<=ROW;i++) {
for(j=1;j<=2;j++) {
U[i][j]=A[i][j];
printf(" %+lf",U[i][j]);
}
printf("\n");
}
for(i=1;i<=ROW-1;i++) {
for(j=ROW-i,t=1;j>=1;j--,t++) {
L[t][i]=U[t+1][2]-U[t][2];
}
for(t=1;t<=ROW-i;t++)
U[t][2]=L[t][i];
}
differenceTable(A,L,ROW);
printf("\n Value at %+lf by Newton's Forward Interpolation formula is %+lf\n",x,result(A,L,x,ROW));
return 0;
}
double result(double **A,double **L,double x,int ROW) {
int i,j;
double h,value=0,p,tmp=1.0;
h=A[2][1]-A[1][1];
p=(x-A[1][1])/h;
for(i=0;i<ROW;i++){
if(i==0)
value +=A[1][2];
else {
for(j=1;j<=i;j++)
tmp=tmp*(p-j+1);
tmp = tmp/factorial(i);
tmp = tmp*L[1][i];
value +=tmp;
}
if(i==0)
printf("%2.2lf ",A[1][2]);
else
printf("%2.2lf ",L[1][i]);
tmp=1.0;
}
return value;
}
double factorial(int n) {
if(n<=1)
return 1;
else
return n*factorial(n-1);
}
void differenceTable(double **A,double **L,int ROW) {
int **position,*ap,j,m,i;
position=(int **)malloc((ROW+1)*sizeof(int*));
ap=(int *)malloc((ROW+1)*sizeof(int));
for(i=0;i<=ROW;i++)
position[i]=(int *)malloc((ROW+1)*sizeof(int));
int an,tmp;
tmp=ROW;
for(i=1;i<=ROW;i++) {
an=1+(i-1)*2;
ap[i]=an;
for(j=1;j<=tmp;j++) {
position[i][j]=an+(j-1)*4;
}
tmp--;
}
tmp=ROW;
tmp=ROW+(ROW-1)*3;
int *match,tmp1,l,*pos,flag,z,k;
match=(int *)malloc((ROW+1)*sizeof(int));
pos=(int *)malloc((ROW+1)*sizeof(int));
for(i=0;i<=ROW;i++) {
match[i]=0;
pos[i]=0;
}
printf("\n-----------------------------Difference Table-----------------------------\n");
printf("\n X Y");
for(i=1;i<=ROW-1;i++)
printf(" D^%d",i);
printf("\n--------------------------------------------------------------------------\n");
for(i=1;i<=tmp;i++) {
tmp1=ROW;
for(l=1;l<=ROW;l++) {
flag=0;
for(m=1;m<=tmp1;m++) {
if(i==position[l][m]) {
flag=1;
match[l]=1;
pos[m]=position[l][m];
break;
}
} //inner for
if(flag==1) {
for(k=1;k<=ROW;k++) {
if(match[k]==0) {
printf("");
}
else {
if(k==1)
//printf("(%d,%d)+(%d,%d)|%d",(position[l][m]/4)+1,k,(position[l][m]/4)+1,k+1,i);
printf("%+2.2lf\t%+2.2lf",A[(position[l][m]/4)+1][k],A[(position[l][m]/4)+1][k+1]);
else {
z=(position[l][m]-ap[k])/4+1;
//printf("\t\t(%d,%d)|%d",z,k-1,i);
printf("\t\t%+2.2lf",L[z][k-1]);
}
}
} //end of k-loop
} //end of flag==1
else
{printf(" ");}
tmp1--;
for(k=0;k<=ROW;k++) {
pos[k]=0;
match[k]=0;
}
}
printf("\n");
}
printf("\n--------------------------------------------------------------------------\n");
}
No comments:
Post a Comment