shuzhi2.cpp

采用线性插值

#include "stdio.h"
#include "math.h"

void main()
{
	double x0[14]={0.00,0.00,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00};/*原始数据x*/
	double y[14]={0.00,0.00,5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00};/*原始数据y*/
	double xx1[200];                                                          /*定义数组x 用此存放按单位变化的x*/
	double yy1[200],yy2[200],yy3[200],yy4[200];                             /*定义数组y,用来存放相应的y*/
	double x,px,middle,h0,px1,px2,px3,px4,hi,hi1;
   	int i,j,k,n;
	double a[13][13];                              /*用于Doolittle分解法求解方程组*/
    double alfa[13],gma[13],beta[13],l[13],m0[13],p[13],q[13],r[13],s[13];                         /*用于三次样条插值法*/
	
	/*********************************************分段线性插值**************************************************/
	  for (i=1;i<192;i++)
	  { xx1[i]=i-1;}
	  
	  k=1;
	   for (x=0;x<191;x++)
		{
	     for (i=1;i<14;i++)                        /*确定x 属于哪一段划分区间*/
		 { if (x<x0[i])
		 {n=i-1;
		 break;}
		 }

		px=((x-x0[n+1])*y[n]/(x0[n]-x0[n+1]))+((x-x0[n])*y[n+1]/(x0[n+1]-x0[n]));  /*求出相应区间对应的一次插值多项式 */

        yy1[k]=px;                                          /*得到y值*/
		k++;
		}
	   yy1[k]=0.00;
	/*for (i=1;i<192;i++)
	   {/*printf("%.2f,%.2f  ",xx1[i],yy1[i]);
      
        printf("%.2f  ",yy1[i]);
	  if (i%10==0) printf("\n");
	 
	 }*/


/*********************************分段二次多项式插值*************************************/
      
	 /*for (i=2;i<14;i++)
	  {x1[i-1]=(x0[i-1]+x0[i])/2;
	   y1[i-1]=(y[i-1]+y[i])/2;}
	  
	   k=1;
	   for (x=0;x<191;x++)
		{
	     for (i=2;i<12;i++)
		 { 
			 if (x<=x1[i])
		      {n=i-1;
		      break;}
			 else n=11;
		 }
        px1=((x-x0[n+1])*(x-x1[n+1])*y1[n]/((x1[n]-x0[n+1])*(x1[n]-x1[n+1])));
		px2=((x-x1[n])*(x-x1[n+1])*y[n+1]/((x0[n+1]-x1[n])*(x0[n+1]-x1[n+1])));
		px3=((x-x1[n])*(x-x0[n+1])*y1[n+1]/((x1[n+1]-x1[n])*(x1[n+1]-x0[n+1])));
		px=px1+px2+px3;
        yy2[k]=px;
		k++;
		}
	/*  for (i=1;i<192;i++)
	   {printf("%.2f,%.2f  ",xx1[i],yy2[i]);}
	   */

       k=1;
	   for (x=0;x<191;x++)
		{
	     for (i=2;i<14;i++)
		 {
		 if (x<=(x0[1]+x0[2])/2) n=1; 
	     else if ((x<=(x0[i]+x0[i+1])/2)&(x>(x0[i-1]+x0[i])/2))                 /*确定x 属于哪一段划分区间*/
		          {n=i-1;
		         break;}
		 else n=11;
		 }
        px1=(x-x0[n+1])*(x-x0[n+2])*y[n]/((x0[n]-x0[n+1])*(x0[n]-x0[n+2]));
		px2=(x-x0[n])*(x-x0[n+2])*y[n+1]/((x0[n+1]-x0[n])*(x0[n+1]-x0[n+2]));
		px3=(x-x0[n])*(x-x0[n+1])*y[n+2]/((x0[n+2]-x0[n])*(x0[n+2]-x0[n+1]));
		px=px1+px2+px3;                                                        /*求出相应区间对应的二次插值多项式 */


        yy2[k]=px;                                                          /*得到y值*/
		k++;
		}
		/*for (i=1;i<192;i++)
	   {/*printf("(%.2f,%.2f)  ",xx1[i],yy2[i]);
	    printf("%.2f  ",yy2[i]);}*/


	   
/********************************分段三次插值**********************************************/

 k=1;
	   for (x=0;x<191;x++)
		{
	     for (i=3;i<13;i++)
		 { if (x<x0[i])
		 {n=i-2;
		 break;}
		 else n=10;                              /*确定x 属于哪一段划分区间*/
		 }
        px1=(x-x0[n+1])*(x-x0[n+2])*(x-x0[n+3])*y[n]/((x0[n]-x0[n+1])*(x0[n]-x0[n+2])*(x0[n]-x0[n+3]));
		px2=(x-x0[n])*(x-x0[n+2])*(x-x0[n+3])*y[n+1]/((x0[n+1]-x0[n])*(x0[n+1]-x0[n+2])*(x0[n+1]-x0[n+3]));
		px3=(x-x0[n])*(x-x0[n+1])*(x-x0[n+3])*y[n+2]/((x0[n+2]-x0[n])*(x0[n+2]-x0[n+1])*(x0[n+2]-x0[n+3]));
		px4=(x-x0[n])*(x-x0[n+1])*(x-x0[n+2])*y[n+3]/((x0[n+3]-x0[n])*(x0[n+3]-x0[n+1])*(x0[n+3]-x0[n+2]));
		px=px1+px2+px3+px4;                         /*求出相应区间对应的二次插值多项式 */

        yy3[k]=px;                                       /*得到y值*/
		k++;
		}
	  	/*for (i=1;i<192;i++)
	   {printf("%.2f  ",yy3[i]);}*/


/******************************三次样条插值************************************************/
for (i=0;i<13;i++)
{
	for (j=0;j<13;j++)
	{a[i][j]=0.0;}                   /*初始化系数矩阵*/
}

for (i=1;i<12;i++)
{hi=x0[i+1]-x0[i];
 hi1=x0[i+2]-x0[i+1];
alfa[i]=hi1/(hi+hi1);
gma[i]=1-alfa[i];
beta[i]=(6/(hi+hi1))*(((y[i+2]-y[i+1])/hi1)-(y[i+1]-y[i])/hi);  /*求出矩阵中的alfa,beta,gma值*/

}

alfa[12]=(x0[2]-x0[1])/(x0[2]-x0[1]+x0[13]-x0[12]);
gma[12]=1-alfa[12];
beta[12]=(6/(x0[2]-x0[1]+x0[13]-x0[12]))*((y[2]-y[1])/(x0[2]-x0[1])-(y[13]-y[12])/(x0[13]-x0[12]));


for (i=1;i<13;i++)
{
	a[i][i+1]=alfa[i];
    a[i][i-1]=gma[i];
	for (j=1;j<13;j++)
	{
		
		if (i==j) {a[i][j]=2;}                           /*得到系数矩阵和右端项*/
		}
	}
	a[1][12]=gma[1];
	a[12][1]=alfa[12];
	
  /*for (i=1;i<13;i++)
  {
	  for (j=1;j<13;j++)
	  {printf("%f  ",a[i][j]);}
	  printf("\n");
  }
  for (i=1;i<13;i++)
  {printf("%f  ",beta[i]);}
  printf("\n");*/
/*********************利用Doolittle分解法反解方程求得新的u0********************************************/
/******分解过程***********/
/*for (j=1;j<13;j++)
{
	u1[1][j]=a[1][j];
}

for (i=2;i<13;i++)
{
	l1[i][1]=a[i][1]/u1[1][1];
}


for (im=1;im<13;im++)
{
	for (j=im;j<13;j++)
	{

		middle=0.0;
		for (tt=1;tt<im;tt++)
		{ middle=middle+l1[im][tt]*u1[tt][j];}
		u1[im][j]=a[im][j]-middle;
	}

	if (im<12)
	{
		for (i=im+1;i<13;i++)
		{
			middle=0.0;
			for (tt=1;tt<im;tt++)
			{
		    middle=middle+l1[i][tt]*u1[tt][im];}
			l1[i][im]=(a[i][im]-middle)/u1[im][im];
		}

	}
}*/



/*************分解过程结束************************************************/
/*****求解过程*****/
/*l[1]=beta[1];
for (i=2;i<13;i++)
{
	middle=0.0;
	for (tt=1;tt<i;tt++)
	{middle=middle+l1[i][tt]*l[tt];}

	l[i]=beta[i]-middle;}

m0[12]=l[12]/u1[12][12];

for (i=11;i>0;i--)
{
	middle=0.0;
	for (tt=i+1;tt<13;tt++)
	{middle=middle+u1[i][tt]*m0[tt];}
	m0[i]=(l[i]-middle)/u1[i][i];
}
m0[0]=m0[12];
for (i=0;i<13;i++)
{printf("%f  ",m0[i]);
}
printf("\n");*/



n=12;
p[1]=2;
for (i=1;i<n-1;i++)
{ q[i]=alfa[i]/p[i];
p[i+1]=2-gma[i+1]*q[i];}

s[1]=gma[1]/p[1];
for (i=2;i<n-1;i++)
{ s[i]=-gma[i]*s[i-1]/p[i];
}
s[11]=(alfa[11]-gma[11]*s[10])/p[11];

r[1]=alfa[12];
for (j=2;j<n-1;j++)
{ r[j]=-r[j-1]*q[j-1];}
r[n-1]=gma[n-1]-r[n-2]*q[n-2];
middle=0.0;
for (j=1;j<n;j++)
{ middle=middle+r[j]*s[j];}
r[n]=2-middle;
  

l[1]=beta[1]/p[1];
for (i=2;i<n;i++)
{ l[i]=(beta[i]-gma[i]*l[i-1])/p[i];
}
middle=0.0;
for (j=1;j<n;j++)
{ middle=middle+r[j]*l[j];}
l[n]=(beta[n]-middle)/r[n];

m0[n]=l[n];
m0[n-1]=l[n-1]-s[n-1]*m0[n];
for (i=n-2;i>0;i--)
{ m0[i]=l[i]-q[i]*m0[i+1]-s[i]*m0[n];}

m0[0]=m0[12];
for (i=0;i<13;i++)
{printf("%f  ",m0[i]);
}
printf("\n");
/****求解过程结束，得到m0************/

       k=1;
	   for (x=0;x<191;x++)
		{
	         for (i=1;i<14;i++)
		 { 
			 if (x<x0[i])
		      {n=i-1;
		      break;}
			 		 }
		h0=x0[n+1]-x0[n];

		/*px=m0[n]*(x0[n+1]-x)*(x0[n+1]-x)*(x0[n+1]-x)/(6*h0)+m0[n+1]*(x-x0[n])*(x-x0[n])*(x-x0[n])/(6*h0)+(y[n]/h0-m0[n]*h0/6)*(x0[n+1]-x)+(y[n+1]/h0-m0[n+1]*h0/6)*(x-x0[n]);/*得到s*/
        px=m0[n-1]*(x0[n+1]-x)*(x0[n+1]-x)*(x0[n+1]-x)/(6*h0)+m0[n]*(x-x0[n])*(x-x0[n])*(x-x0[n])/(6*h0)+(y[n]/h0-m0[n-1]*h0/6)*(x0[n+1]-x)+(y[n+1]/h0-m0[n]*h0/6)*(x-x0[n]);
        yy4[k]=px;/*得到y*/
		k++;
		}
for (i=1;i<192;i++)
	   {printf("%.2f  ",yy4[i]);
if (i%10==0) printf("\n");}
	   

}
/*******************************主程序结束*******************/