nihe.cpp

对给定数据进行最小二乘曲面拟和

#include<stdio.h>
#include<math.h>
#include <malloc.h>
#define UT  6
#define MAXK 11
#define POINTX 11//x点数
#define POINTY 21//y点数
#define MM    30  
FILE *fp;



int  max(double s[],int k,int n)
{ int q,i;
	q=k;
  for (i=k+1;i<=n;i++)
   if(fabs(s[i])>fabs(s[k]))
     q=i;
   return(q);
printf("\n %d  ",q);
}

void change(double  low[][MM+1],double  a[][MM+1],double s[],double b[],int k,int q,int n)
{int t;
 double  temp,temp1;
  if(q!=k)
	  { if(k==1)
        {temp1=s[k];
         s[k]=s[q];
         s[q]=temp1;
         }
	  
	  for(t=1;t<=k-1;t++)
     {temp=low[k][t];
      low[k][t]=low[q][t];
      low[q][t]=temp;
      temp1=s[k];
      s[k]=s[q];
      s[q]=temp1;
       }
     for (t=k;t<=n;t++)
     {temp=a[k][t];
      a[k][t]=a[q][t];
      a[q][t]=temp;
      }
     if(k<n)
      {temp=b[k];
      b[k]=b[q];
      b[q]=temp;
       }
     }
}
void count(double low[][MM+1],double a[][MM+1],double u[][MM+1],double s[],int k,int n)
{int j,t,i;
  double sum;
u[k][k]=s[k];
 for(j=k+1;j<=n;j++)
   {sum=0;
     if(k>1)
     for (t=1;t<=k-1;t++)
      sum=sum+low[k][t]*u[t][j];
     
     u[k][j]=a[k][j]-sum;
     for(i=k+1;i<=n;i++)
     low[i][k]=s[i]/u[k][k];
    }
}

void ans(double y[],double b[],double x[],double low[][MM+1],double u[][MM+1],int k,int n)
{ int i,t;
  double sum;
  y[1]=b[1];
 for(i=2;i<=n;i++)
   {sum=0;
    for(t=1;t<=i-1;t++)
     sum=sum+low[i][t]*y[t];
   y[i]=b[i]-sum;
    }
  x[n]=y[n]/u[n][n];
  printf("\nou:%lf",x[n]);
  for(i=n-1;i>=1;i--)
   {sum=0;
    for(t=i+1;t<=n;t++)
     sum=sum+u[i][t]*x[t];
    x[i]=(y[i]-sum)/u[i][i];
   }
}

 
void change_dl(double a[][M+1],double b[],double x[],int n)

{ int i,j,k,q,t,or,qu,n;     
  double sum,low[MM+1][MM+1],u[MM+1][MM+1],s[MM+1],y[MM+1];
  int max(double s[],int k,int n);
void enter(double  a[][MM+1],double  b[],int n);
void change(double  low[][MM+1],double  a[][MM+1],double s[],double b[],int k,int q,int n);
  void count(double low[][MM+1],double a[][MM+1],double u[][MM+1],double s[],int k,int n);
  void ans(double y[],double b[],double x[],double low[][MM+1],double u[][MM+1],int k,int n);
//n=M;
//FILE *fp;
  //if((fp=fopen("E:\luyan.txt","w"))==NULL)
  // {printf("can't open file");
    //}

//printf("\nenter n:");
//scanf("%d",&n);
//enter(a,b,n);
for(i=1;i<=n;i++)
  for(j=1;j<=n;j++)
   {u[1][j]=a[j][1];
  low[i][1]=a[i][1]/u[1][1];
    }
for (k=1;k<=n;k++)
  { for (i=k;i<=n;i++)
      {  sum=0;
      if(k>1)
        for (t=1;t<=k-1;t++)
       sum=sum+low[i][t]*u[t][k];
		s[i]=a[i][k]-sum;
     }
      q=max(s,k,n);
      printf("\n %d  ",q);
       for (or=k;or<=n;or++)
		  printf("\ns[%d]: %lf  ",or,s[or]);
      
	  
      change(low,a,s,b,k,q,n);
      for (or=k;or<=n;or++)
		  printf("\ns[%d]: %lf  ",or,s[or]);
	  

      count(low,a,u,s,k,n);

	  for (qu=k+1;qu<=n;qu++)
        printf("\nu: %lf  low:%lf ",u[k][qu],low[qu][k]);
  }
 ans(y,b,x,low,u,k,n);
printf("\n");
 for (i=1;i<=n;i++)
 {//fprintf(fp,"x[%d]: %lf  ",i,x[i]);
   printf("y[%d]: %lf  \n",i,y[i]);
}
 printf("\n u:\n");
for(i=1;i<=n;i++)
{for(j=1;j<=n;j++)
{ if(i>j) u[i][j]=0;
   printf("%lf    ",u[i][j]);
}
  printf("\n");
}
 printf("\n l:\n");
for(i=1;i<=n;i++)
{for(j=1;j<=n;j++)
{ if(j>=i) low[i][j]=0;
   printf("%lf    ",low[i][j]);
    //fprintf(fp,"%lf    ",low[i][j]);
	//if(j%3==0)
    //fprintf(fp,"\n");
}
  printf("\n");
}


}





//曲面拟和
   pstar(int k)
   {

	   double matb[MM+1][MM+1],matalfa[MM+1],matg[MM+1][MM+1],matbtb[MM+1][MM+1];
       double matgtg[MM+1][MM+1];
      for(i=1;i<=POINTX;i++)
		  for(j=1;j<=k+1;j++)
			  matb[i][j]=pow(x[i],j-1);//求B矩阵
      
      for(i=1;i<=POINTY;i++)
		  for(j=1;j<=k+1;j++)
			  matg[i][j]=pow(y[i],j-1);//求G矩阵
		  //求BT*B
      for(i=1;i<=k+1;i++)
		  for(j=1;j<=k+1;j++)
		  {   matbtb[i][j]=0;
			  for(mm=1;mm<=POINTX;mm++)
			   matbtb[i][j]=matb[mm][i]*matb[mm][j]+matbtb[i][j];
		  }
	  //求GT*G
      for(i=1;i<=k+1;i++)
		  for(j=1;j<=k+1;j++)
		  {   matgtg[i][j]=0;
			  for(mm=1;mm<=POINTX;mm++)
			   matgtg[i][j]=matg[mm][i]*matg[mm][j]+matgtg[i][j];
		  }
        //计算BT*uj
        for(j=1;j<=POINTY;j++)
		{
			for(i=1;i<=k+1;i++)
		  {   uj[i]=0;
			 for(mm=1;mm<=POINTX;m++)
			  uj[i]=matb[mm][i]*u[mm][1]+uj[i];
			}   //计算BT*uj
          change_dl(matbtb,uj,matalfa,k+1);
          

             qj[j]=0;
          for(i=1;i<=k+1;i++)
		   qj[j]=qj[j]+matalfa[i]*pow(x0,i-1);//x0为给定插值点
		}










 void main()
	{ int n;
      double t0,u0,fz0;
     void tqip( int n1,int m1,double *x,double *y,double *z,double u,double v, double *f );
	 dataread();
     t0=0.243185;
	 u0=1.345231;
     tqip(UT,UT,&t[1],&u[1],&fz[1],t0,u0,&fz0);
      printf("fz0=%lf\n",fz0);
      scanf("%d",&n);
	}