analysis2.txt

数值分析B计算实习作业二：分别用分段线性插值、分段二次多项式插值、 分段三次多项式插值和三次样条插值对所给的数据进行细化

/*------------------------------------------
数值分析B计算实习作业二：分别用分段线性插值、分段二次多项式插值、
分段三次多项式插值和三次样条插值对所给的数据进行细化

-------------------------------------------*/
#include<stdio.h>
#include<math.h>
#define N 12
double x[13]={0.0,4.74,9.5,19.0,38.0,57.0,76.0,95.0,114.0,133.0,152.0,171.0,190.0}; 
double y[13]={0.0,5.32,8.10,11.97,16.15,17.1,16.34,14.63,12.16,9.69,7.03,3.99,0.00}; //机翼轮廓线部分数据
double M[13];  //三次样条插值法中，各插值节点的二次导数值
double or_va[4];
/*------------------------------------------
子函数功能：分段线性插值
入口参数：区间[X0,Xn]之间的整数m
输出参数：分段线性插值法细分后对应于m的函数值
--------------------------------------------*/
double order_1(int m)  //order_1() 分段线性插值
{
	int i;
	double k,v;
	if(m>=0.0&&m<=190.0)
	{
		for(i=0;i<12;i++)
		{
		    k=(y[i+1]-y[i])/(x[i+1]-x[i]);
		    if(m>=x[i]&&m<=x[i+1])
			{
                v=y[i]+k*(m-x[i]);  //在区间[Xi，Xi+1]间一次多项式
			    return(v);
			}
		    else continue;
		}
	}
	else return(-1);
}
/*------------------------------------------
子函数功能：分段二次多项式插值
入口参数：区间[X0,Xn]之间的整数m
输出参数：分段二次多项式插值法细分后对应于m的函数值
--------------------------------------------*/
double order_2(int m)  
{
	int i;
	double a[12],b[12],c[12],v,dit,dita,ditb,ditc,d1,d2,d3;
    for(i=1;i<12;i++)
		{
		    d1=pow(x[i-1],2);
            d2=pow(x[i],2);
			d3=pow(x[i+1],2);
		    dit=(x[i]-x[i-1])*(x[i+1]-x[i-1])*(x[i+1]-x[i]);
            dita=y[i-1]*(x[i]-x[i+1])-y[i]*(x[i-1]-x[i+1])+y[i+1]*(x[i-1]-x[i]);
			ditb=y[i]*(d1-d3)-y[i-1]*(d2-d3)-y[i+1]*(d1-d2);
			ditc=y[i-1]*(d2*x[i+1]-d3*x[i])-y[i]*(d1*x[i+1]-d3*x[i-1])+y[i+1]*(d1*x[i]-d2*x[i-1]);
			a[i]=-dita/dit;
			b[i]=-ditb/dit;
			c[i]=-ditc/dit;
		    
		}
	while(m>=x[0]&&m<(x[1]+x[0])/2) 
	{v=a[1]*pow(m,2)+b[1]*m+c[1];
			    return(v);}
	while(m>=(x[11]+x[12])/2&&m<=x[12])
	{v=a[11]*pow(m,2)+b[11]*m+c[11];
			    return(v);}
	for(i=1;i<12;i++)
	{
         if(m>=(x[i]+x[i-1])/2&&m<(x[i+1]+x[i])/2)
			{
                v=a[i]*pow(m,2)+b[i]*m+c[i];
			    return(v);
			}
		    else continue;
	}
}
/*-------------------------------------------
子程序功能：高斯消元法求线性方程组Ax=b的解
入口参数：系数矩阵a[][]，b[]；
无输出参数
-------------------------------------------*/
void gauss(double a[4][4],double b[4])
{
	int i,j,k;
	double mik,orvass;
	for(k=0;k<3;k++)
	{
		for(i=k+1;i<4;i++)
		{
			mik=a[i][k]/a[k][k];
			for(j=0;j<4;j++)
			{a[i][j]-=mik*a[k][j];}
            b[i]-=mik*b[k];
		}
	}
	or_va[3]=b[3]/a[3][3];
	for(k=2;k>=0;k--)
	{
		orvass=0;
		for(j=k+1;j<4;j++)
		{orvass+=a[k][j]*or_va[j];}
		or_va[k]=(b[k]-orvass)/a[k][k];
	}
}
/*------------------------------------------
子函数功能：分段三次多项式插值
入口参数：区间[X0,Xn]之间的整数m
输出参数：分段三次多项式插值法细分后对应于m的函数值
--------------------------------------------*/
double order_3(int m)
{
	int i,j,k;
	double v,a[10],b[10],c[10],d[10],A[4][4],be[4];
    for(k=1;k<11;k++)
	{
        for(i=0;i<4;i++)
		{
			for(j=0;j<4;j++)
			{
				A[i][j]=pow(x[k+2-i],3-j);
			}
			be[i]=y[k+2-i];
		}
		gauss(A,be);
		a[k]=or_va[0];
		b[k]=or_va[1];
		c[k]=or_va[2];
		d[k]=or_va[3];
	}
while(m>=x[0]&&m<x[2]) 
	{v=a[1]*pow(m,3)+b[1]*pow(m,2)+c[1]*m+d[1];
			    return(v);}
	while(m>=x[10]&&m<=x[12])
	{v=a[10]*pow(m,3)+b[10]*pow(m,2)+c[10]*m+d[10];
			    return(v);}
	for(i=2;i<10;i++)
	{
         if(m>=x[i]&&m<x[i+1])
			{
                v=a[i]*pow(m,3)+b[i]*pow(m,2)+c[i]*m+d[i];
			    return(v);
			}
		    else continue;
	}
}
/*------------------------------------------
子函数功能：追赶法求拟三对角线性方程组的解
入口参数：组成拟上对角矩阵的向量a[N],b[N],c[N],d[N],组成参见书本2.2.5节
--------------------------------------------*/
void pursue(double a[N],double c[N],double d[N],double b[N])
{
	int i;
	double rss,yss,pur_p[N-1],pur_q[N-1],pur_s[N-1],pur_r[N],pur_y[N];
	pur_p[0]=a[0];
	for(i=0;i<N-2;i++)
	{
		pur_q[i]=c[i]/pur_p[i];
		pur_p[i+1]=a[i+1]-d[i+1]*pur_q[i];
	}
	pur_s[0]=d[0]/pur_p[0];
    for(i=1;i<N-2;i++)
	{
		pur_s[i]=-d[i]*pur_s[i-1]/pur_p[i];
	}
	pur_s[N-2]=(c[N-2]-d[N-2]*pur_s[N-3])/pur_p[N-2];
	pur_r[0]=c[N-1];
    for(i=1;i<N-2;i++)
	{
		pur_r[i]=-pur_q[i-1]*pur_r[i-1];
	}
	pur_r[N-2]=d[N-1]-pur_r[N-3]*pur_q[N-3];
	rss=0;
	for(i=0;i<N-1;i++)
	{
		rss+=pur_r[i]*pur_s[i];
	}
	pur_r[N-1]=a[N-1]-rss;
	pur_y[0]=b[0]/pur_p[0];
    for(i=1;i<N-1;i++)
	{
		pur_y[i]=(b[i]-d[i]*pur_y[i-1])/pur_p[i];
	}
	yss=0;
	for(i=0;i<N-1;i++)
	{
		yss+=pur_r[i]*pur_y[i];
	}
	pur_y[N-1]=(b[N-1]-yss)/pur_r[N-1];
	M[N-1]=pur_y[N-1];
	M[N-2]=pur_y[N-2]-pur_s[N-2]*M[N-1];
	for(i=N-3;i>=0;i--)
	{
		M[i]=pur_y[i]-pur_q[i]*M[i+1]-pur_s[i]*M[N-1];
	}
}
/*------------------------------------------
子函数功能：三次样条插值
入口参数：区间[X0,Xn]之间的整数m
输出参数：三次样条插值法细分后对应于m的函数值
--------------------------------------------*/
double order_4(int n)
{
	int i;
	double v,aefa[12],beta[12],gama[12],cont[12],k1,k2,h1,h2;
	for(i=0;i<11;i++)
	{
		h1=x[i+1]-x[i];h2=x[i+2]-x[i+1];
		aefa[i]=h2/(h1+h2);
		beta[i]=6*((y[i+2]-y[i+1])/h2-(y[i+1]-y[i])/h1)/(h1+h2);
	}
	aefa[11]=(x[1]-x[0])/(x[1]-x[0]+x[12]-x[11]);
	beta[11]=6*((y[1]-y[0])/(x[1]-x[0])-(y[12]-y[11])/(x[12]-x[11]))/(x[12]-x[11]+x[1]-x[0]);
	for(i=0;i<12;i++)
	{
		cont[i]=2;
		gama[i]=1-aefa[i];
	}
    pursue(cont,aefa,gama,beta);
	for(i=12;i>0;i--)
	{M[i]=M[i-1];}
	M[0]=M[11];
		for(i=0;i<12;i++)
		{
		    if(n>=x[i]&&n<=x[i+1])
			{
                h1=x[i+1]-x[i];
                k1=y[i]/h1-(M[i]*h1)/6;
		        k2=y[i+1]/h1-(M[i+1]*h1)/6;
                v=(M[i]*pow(x[i+1]-n,3)+M[i+1]*pow(n-x[i],3))/(6*h1)+k1*(x[i+1]-n)+k2*(n-x[i]);
			    return(v);
			}
		    else continue;
		}
}
/*-----------------------------------
主函数功能：调用相应的子程序对数据用分段线性、分段二次多项式、
分段三次多项式和三次样条插值法进行细分并输出相应的结果到yyj.txt
--------------------------------------*/
void main()
{
	int i;
	double result_1[191],result_2[191],result_3[191],result_4[191];
	FILE *fish;
	fish=fopen("D:\\yyj.txt","w");
	for(i=0;i<=190;i++)
	{
         result_1[i]=order_1(i);
		 result_2[i]=order_2(i);
		 result_3[i]=order_3(i);
		 result_4[i]=order_4(i);
		 fprintf(fish,"%d   %f   %f   %f   %f\n",i,result_1[i],result_2[i],result_3[i],result_4[i]);
	}
	fclose(fish);
	
}