Pages

Tuesday, 11 October 2011

Gauss Backward Difference Formula

#include<stdio.h>
#include<math.h>
float factorial(int);
void main()
{
        float ax[10],y[10];
        float dy[10][10],u,h,x,a,sum=0,l,oddprod=1,evenprod=1;
        int i,j,k,n,mid=0;
        printf("\n\t\tGauss Forward Difference Formula");
        printf("\n=======================================================\n");
        printf("Enter the number of terms ");
        scanf("%d",&n);
        printf("\n=======================================================\n");
        printf("Enter the array ax\n");
        for(i=0;i<n;i++)
        {
                printf("Enter the value of ax[%d] = ",i);
                scanf("%f",&ax[i]);
        }
        printf("\n=======================================================\n");
        for(i=0;i<n;i++)
        {
                printf("Enter the value of y[%d] = ",i);
                scanf("%f",&y[i]);
                
        }
        printf("\nEnter the interpolation point ");
        scanf("%f",&x);
        h=ax[1]-ax[0];
        for(i=0;i<n;i++)
        {
                for(j=0;j<=n;j++)
                        dy[i][j]=-1;
                
        }
        for(i=0;i<n;i++)
        {
                dy[i][0]=y[i];
        }
        
        k=1;
        j=1;
        while(j<n)
        {
                i=1;
                while(i<=n-k)
                {
                        dy[i-1][j]=dy[i][j-1]-dy[i-1][j-1];
                        i++;
                }
                j++;
                k++;
        }
        i=0;
        printf("\n================================================================================\n");
        printf("\ty");
        for(i=1;i<n;i++)
                printf("\td^%dy",i);
        printf("\n");
        for(i=0;i<n;i++)
        {
                for(j=0;j<n;j++)
                {
                        if(dy[i][j]==-1.0)
                                printf("\t_");
                        else
                                printf("\t%-6.3f",dy[i][j]);
                }
                printf("\n");
        }
        i=0;
        do{
                i++;
        }while(ax[i]<x);
        mid=i;
        u=(x-ax[mid])/h;
        sum=sum+y[mid]+(u*dy[mid-1][1]);
        mid=mid-1;
        for(i=2;i<n;i++)
        {
                
                if(i%2==1)
                {
                        oddprod=1;
                        for(l=1;l<=i/2;l++)
                        {
                                oddprod=oddprod*(u-l);
                                        
                        }
                        mid-=1;
                        sum=sum + (((evenprod*u*oddprod)/factorial(i))*dy[mid][i]);
                        
                }
                else
                {
                        evenprod=1;
                        for(l=1;l<=i/2;l++)
                        {
                                evenprod=evenprod*(u+l);
                                        
                        }
                        
                        sum=sum + (((oddprod*u*evenprod)/factorial(i))*dy[mid][i]);
                        
                        
                }
                
        }
        
        printf("\nValue at %6.3f is %6.3f ",x,sum);
        getchar();
}

float factorial(int x)
{
        int i,j=1;
        for(i=1;i<=x;i++)
                j=j*i;
        return j;

}

No comments:

Post a Comment