2order.cpp

二次拟合 线性方程组的求解的函数

  #include   "stdafx.h"   
  #include   "math.h"   
  #include   "stdio.h"   
  #define   ABS(x)   (x)>0?(x):-(x)   
  #define   SWAP(a,b)   {temp=(a);(a)=(b);(b)=temp;}   
    
  void   solve(double   **a,double   *b,double   *x,int   n)   
  {   
  int   i,j,k,ik;   
  double   mik,temp;   
  for(k=0;k<n;k++)   
  {   
  mik=-1;   
  for(i=k;i<n;i++)   
  if(ABS(a[i][k])>mik)   
  {   
  mik=ABS(a[i][k]);   
  ik=i;   
  }   
  for(j=k;j<n;j++)   
  SWAP(a[ik][j],a[k][j]);   
  SWAP(b[k],b[ik]);   
  b[k]/=a[k][k];   
  for(i=n-1;i>=k;i--)   
  a[k][i]/=a[k][k];   
  for(i=k+1;i<n;i++)   
  {   
  b[i]-=a[i][k]*b[k];   
  for(j=n-1;j>=k;j--)   
  a[i][j]-=a[i][k]*a[k][j];   
  }   
  }   
  for(i=n-1;i>=0;i--)   
  {   
  x[i]=b[i];   
  for(j=i+1;j<n;j++)   
  x[i]-=a[i][j]*x[j];   
  }   
  }   
    
  void   linear(double   **x,double   *y,double   *beta,int   n,int   p)   
  {   
  double   **a,*b;   
  int   i,j,k;   
  a=new   double*[n];   
  for(i=0;i<n;i++)   
  a[i]=new   double[p];   
  for(i=0;i<n;i++)   
  for(j=0;j<p;j++)   
  {   
  a[i][j]=0;   
  for(k=0;k<n;k++)   
  a[i][j]+=x[k][i]*x[k][j];   
  }   
  b=new   double[p];   
  for(i=0;i<p;i++)   
  {   
  b[i]=0;   
  for(j=0;j<n;j++)   
  b[i]+=x[j][i]*y[j];   
  }   
  solve(a,b,beta,p);   
  }   
    
  void   polyfit(double   *x,double   *y,double   *beta,int   n,int   p)   
  {   
  int   i,j;   
  double   **xx;   
  xx=new   double*[n];   
  for(i=0;i<n;i++)   
  {   
  xx[i]=new   double[p];   
  for(j=0;j<p;j++)   
  xx[i][j]=pow(x[i],j);   
  }   
  linear(xx,y,beta,n,p);   
  }   
  void   main()   
  {   
  FILE   *fp;   
  int   i,n;   
  double   *x,*y;   
  fp=fopen("d:\\data.txt","r");   
  fscanf(fp,"%d",&n);   
  x=new   double[n];   
  y=new   double[n];   
  for(i=0;i<n;i++)   
  {   
  fscanf(fp,"%lf",x+i);   
  fscanf(fp,"%lf",y+i);   
  }   
  fclose(fp);   
    
  int   p=4;   
  double   *beta;   
  beta=new   double[p];   
  polyfit(x,y,beta,n,p);   
  for(i=0;i<p;i++)   
  printf("%6.4f     ",beta[i]);   
  }