curvefit.c

c语言版的多项式曲线拟合。 用最小二乘法进行曲线拟合. 用p－1 次多项式进行拟合

/*
用最小二乘法进行曲线拟合.
用p－1 次多项式进行拟合，p<= 10;
x,y 的第0个域x[0],y[0],没有用，有效数据从x[1],y[1] 开始
nNodeNum,有效数据节点的个数。
b,为输出的多项式系数，b[i] 为b[i-1]次项。b[0],没有用。
b,有10个元素ok。
*/
void CurveFit(double *x,double *y,int p,int nNodeNum, double *b)
{

	int t=nNodeNum;  
	double a[10][10];
	int i=1;   //核心程序
    int j=1;
    int k=1;
    double ss=0.0;
    int r;
    double m;
    double temp;
	memset(a,0,sizeof a);
    for(j=1;j<=p;j++)
    {
      for(i=1;i<=j;i++)
      {
        ss=0.0;
       for( k=1;k<=t;k++)
        {
          ss=ss+pow(x[k],i+j-2);
        }
        a[i][j] = ss;
        a[j][i]= ss;
      }
    }
    for(i=1;i<=p;i++)

   {
      ss=0.0;
      for(k=1;k<=t;k++)
      {
        ss=ss+y[k]*pow(x[k],i-1);
      }
      b[i] = ss;
    }
    for(k=1;k<=p;k++)
    {
      m=fabs(a[k][k]);r=k;
      //选主元
      for(i=k;i<=p;i++)
      {
        if( m > fabs(a[i][k]) )        
		{
          m=fabs(a[i][k]);
          r=i;
        }
      }
      //交换两行
      if(r>k)
      {
        for(j=k;j<=p;j++)
        {
          temp=a[k][j];
          a[k][j]=a[r][j];
          a[r][j]=temp;
        }
        temp=b[k];
        b[k]=b[r];
        b[r]=temp;
      }
      //消元计算
      for(j=k+1;j<=p;j++)
        a[k][j]=a[k][j]/a[k][k];
      b[k]=b[k]/a[k][k];
      for(i=1;i<=p;i++)
      {
        if(i==k) continue;
        for(j=k+1;j<=p;j++)
          a[i][j]=a[i][j]-a[i][k]*a[k][j];

      }

     for(i=1;i<=p;i++)
      {
        if(i==k) continue;
        b[i]=b[i]-a[i][k]*b[k];

      }
	}
}