nihe4.cpp

四次及四次以下多项式拟和程序,可在c或c++环境下运行

 
#include<math.h>
#define D 0x40000
/*此函数为高斯全主元消去法解n元一次方程组，用于下面拟合函数Guess中计算目标曲线函数的各拟合系数*/

int agaus(double a[ ],double b[ ],int n)
{
 int *js,l,k,i,j,is,p,q;
     double d,t;
     js=new int[n];//(int*)malloc(n*sizeof(int));
     l=1;
  for (k=0;k<=n-2;k++)
  { 
d=0.0;
     for (i=k;i<=n-1;i++)
        for (j=k;j<=n-1;j++)
         { 
t=fabs(a[i*n+j]);
             if (t>d) 
{ d=t; js[k]=j; is=i;}
          }
        if (d+1.0==1.0) l=0;
        else
          { if (js[k]!=k)
              for (i=0;i<=n-1;i++)
                { p=i*n+k; q=i*n+js[k];
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
            if (is!=k)
              { for (j=k;j<=n-1;j++)
                  { p=k*n+j; q=is*n+j;
                    t=a[p]; a[p]=a[q]; a[q]=t;
                  }
                t=b[k]; b[k]=b[is]; b[is]=t;
              }
          }
        if (l==0)
          { 
			delete []js; //printf("fail\n");
            return(0);
          }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
          { p=k*n+j; a[p]=a[p]/d;}
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)
          { for (j=k+1;j<=n-1;j++)
              { p=i*n+j;
                a[p]=a[p]-a[i*n+k]*a[k*n+j];
              }
            b[i]=b[i]-a[i*n+k]*b[k];
          }
      }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
      { delete []js; //printf("fail\n");
        return(0);
      }
    b[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
      { t=0.0;
        for (j=i+1;j<=n-1;j++)
          t=t+a[i*n+j]*b[j];
        b[i]=b[i]-t;
      }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
    delete []js;
    return(1);
  }







/*此为4次5项式拟合函数，其中数组x，y分别存放将被拟合的数据点的横坐标和纵坐标，
它们的个数为n，g[0],g[1],g[2],g[3],g[4]分别代表拟合后得到的y=a+b*x+c*x*x+d*x*x*x+e*x*x*x*x中的a,b,c */


void nihe4(double x[],double y[],int n,double g[])
{ 
	int i;
	double X=0.0,X2=0.0,X3=0.0,X4=0.0,X5=0.0,X6=0.0,X7=0.0,X8=0.0;
	double S=0.0,Y=0.0,XY=0.0,X2Y=0.0,X3Y=0.0,X4Y=0.0;
	double h[25]={0};
	for(i=0;i<n;i++)
	{
		X=X+x[i];
		X2=X2+x[i]*x[i];
		X3=X3+pow(x[i],3);
		X4=X4+pow(x[i],4);
		X5=X5+pow(x[i],5);
		X6=X6+pow(x[i],6);
		X7=X7+pow(x[i],7);
		X8=X8+pow(x[i],8);
		S=S+1.0;
		Y=Y+y[i];
		XY=XY+x[i]*y[i];
		X2Y=X2Y+x[i]*x[i]*y[i];
		X3Y=X3Y+pow(x[i],3)*y[i];
		X4Y=X4Y+pow(x[i],4)*y[i];
	}
	h[0]=X4;
	h[1]=X5;
	h[2]=X6;
	h[3]=X7;
	h[4]=X8;
	h[5]=X3;
	h[6]=X4;
	h[7]=X5;
	h[8]=X6;
	h[9]=X7;
	h[10]=X2;
	h[11]=X3;
	h[12]=X4;
	h[13]=X5;
	h[14]=X6;
	h[15]=X;
	h[16]=X2;
	h[17]=X3;
	h[18]=X4;
	h[19]=X5;
	h[20]=S;
	h[21]=X;
	h[22]=X2;
	h[23]=X3;
	h[24]=X4;
	g[0]=X4Y;
	g[1]=X3Y;
	g[2]=X2Y;
	g[3]=XY;
	g[4]=Y;
	if(agaus(h,g,5)!=0)
		return;
}