3次样条插值

	# include<stdio.h>  
	# include<math.h>  
	# include<stdlib.h>
	# define delta 1e-6  
	# define N 100  
	int main()  
	{  
	    int i,j,t,r,n,u,c=0;
	    float p,L,max,s,h,D1,D2;
	    float x[N],M[N],y[N];
	    float a[N][N+1]={0};
	    /*输入插值条件*/ 
	    printf("请输入插值节点的个数：\n");    
	    scanf("%d",&n);
	    printf("请从左到右输入插值节点：\n");
	    for(i=0;i<n;i++)
		{
			scanf("%f",&x[i]);
		} 
        printf("请输入插值节点对应的函数值：\n");
	    for(i=0;i<n;i++)
		{
			scanf("%f",&y[i]);
		} 
        printf("请输入插值初始条件：\n");
        printf("y'(0)="); 
        scanf("%f",&D1);
        printf("y'(n)="); 
        scanf("%f",&D2);      
	    /*产生弯矩方程组增广矩阵系数*/ 
	    a[0][1]=1;
	    a[n-1][n-2]=1;
	    for(i=0;i<n;i++)
	        a[i][i]=2;
        for(i=1;i<n-1;i++)
        {
        	h=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        	a[i][i-1]=h;
        	a[i][i+1]=1-h;
        }
        a[0][n]=6*((y[n-2]-y[n-1])/(x[n-2]-x[n-1])-D2)/(x[n-2]-x[n-1]); 
        a[n-1][n]=6*((y[1]-y[0])/(x[1]-x[0])-D1)/(x[1]-x[0]); 
        for(i=1;i<n-1;i++)
        {
        	a[i][n]=6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1]); 
        }
	    printf("请确认弯矩方程的增广矩阵为：\n");
	    for(i=0;i<n;i++)
	    {
	        for(j=0;j<n+1;j++)
	            printf("%3.2f\t",a[i][j]);
	        printf("\n");
	    }
	    for(j=0;j<n-1;j++)
	    {
	        {
	            max=fabs(a[j][j]);
	            r=j;
	        }
	        for(i=j+1;i<n;i++)
	        if(fabs(a[i][j])>max)
	        {
	            max=a[i][j];
	            r=i;
	        }
	        if(max<delta)
	            printf("矩阵奇异\n");
	        for(t=j;t<n+1;t++)
	        {
	             p=a[j][t];
	             a[j][t]=a[r][t];
	             a[r][t]=p;
	        }
	        for(i=j+1;i<N;i++)
	        {
	            L=a[i][j]/a[j][j];
	            for(t=j;t<n+1;t++)
	                a[i][t]=a[i][t]-L*a[j][t];
	        }
	    }
	    /*输出弯矩方程的计算结果*/ 
	    printf("整理后的增广矩阵为：\n");
	    M[n-1]=a[n-1][n]/a[n-1][n-1];
	    for(i=n-2;i>=0;i--)   
	    {
	        s=a[i][n];  
	        for(j=i+1;j<n;j++)  
	        s=s-a[i][j]*M[j];  
	        M[i]=s/a[i][i];   
	    }   
	    for(u=0;u<n;u++)  
	        for(j=0;j<n+1;j++)  
	        {  
	            printf("%3.2f\t",a[u][j]);  
	            c++;  
	            if(c%(n+1)==0) 
	            printf("\n");  
	        }  
        printf("弯矩方程解向量为：\n"); 
	    for(i=0;i<n;i++)  
	    {    
	        printf("M(%d)=%f\n",i+1,M[i]);  
	        if(i==n-1)  
	        printf("\n");   
	    }  
	    /*计算插值函数*/
	    printf("样条函数如下：\n");
		for(i=0;i<n-1;i++)
		{
			printf("S(x)=%f+%f(x-%.2f)+%f(x-%.2f)^2+%f(x-%.2f)^3\t%2.1f<x<%2.1f",y[i],((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M[i+1]/6)*(x[i+1]-x[i])),x[i],M[i]/2,x[i],(M[i+1]-M[i])/6/(x[i+1]-x[i]),x[i],x[i],x[i+1]);
			printf("\n");
		} 
		printf("S(i+0.5)的值如下：\n");
		for(i=0;i<n-1;i++)
		{
			printf("S(%d+0.5)=%f",i,y[i]+((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M[i+1]/6)*(x[i+1]-x[i]))*(i+0.5-x[i])+M[i]/2*(i+0.5-x[i])*(i+0.5-x[i])+(M[i+1]-M[i])/6/(x[i+1]-x[i])*(i+0.5-x[i])*(i+0.5-x[i])*(i+0.5-x[i]));
			printf("\n");
		} 
	    system("pause");  
	    return 0;  
	}