leastsquarfit.txt

用最小二乘法拟合曲线y=a0+a1*x+a2*x^2+a3*x^3+...+an*x^n 的vc源码

//--------------------------------------------------------------------------- 
//用最小二乘法拟合曲线y=a0+a1*x+a2*x^2+a3*x^3+...+an*x^n; 
//x,y为自变量和因变量数组;n为为系数个数;m为数据的组数;sum用于存和的数组; 
//l为sum数组边界;r为系数数组,也用于存x^n*y的和;a用于存放正则方程组的系数方阵; 
//返回值:0,正常;-1,正则方程系数方阵异常;-2,数据太少,无法拟合; 
int LeastSquarFit(double x[],double y[],int n,int m,double sum[], 
                double an[],double dMatrix[][5]) 
{ 
  int i,j,k,k1,rtn=0; 
  double w,dx,dy,xx,eps; 
        if(m<n) //数据的组数小于所求系数的个数,无法拟合曲线! 
                return -2; 
        sum[0]=m; 
        for(j=0;j<n;j++) 
                an[j]=0.0; 
        for(j=1;j<2*n-1;j++) 
                sum[j]=0.0; 
        for(i=0;i<m;i++) 
        { 
                dx=1.0; 
                dy=y[i]; 
                xx=x[i]; 
                for(j=1;j<2*n-1;j++) 
                { 
                        dx*=xx; 
                        sum[j]+=dx; //x的j次方的和 
                } 
                for(j=0;j<n;j++) 
                { 
                        an[j]+=dy; //y与x的j次方的乘积的和 
                        dy*=xx; 
                } 
        } 
        k1=0; 
        for(j=0;j<n;j++) 
        { 
                k=k1; 
                for(i=0;i<n;i++) 
                { 
                        dMatrix[i][j]=sum[k]; //正则方程的系数为各次x的和 
                        k++; 
                } 
                k1++; 
        } 
        int nans[5]; 
        rtn=SolvLinEqua(dMatrix,n,an,nans,1e-10); 
        return rtn; 
} 
//--------------------------------------------------------------------------- 
//用高斯消除法解n阶一次线性方程组,先化方阵为对角阵,然后求解 
//Matrix为系数方阵；n为方阵的阶数；an作为输入参数时为方程组右端常数项， 
//an作为返回值时为方程组的解；min为误差限，最小为1e-8；nans为是计算中要用的中间数组 
//函数返回值：0为正常；-1，系数矩阵异常 
int SolvLinEqua(double dMatrix[][5],int n,double an[],int nans[],double min) 
{ 
        int i,k,j,m,nstop; 
        double epss,bigst,pivot,v,w; 
        nstop=0; //返回值预置为0-正常 
        if(min<1e-8) //误差限至少为1e-8 
                min=1e-5; 
        epss=min*min; 
        for(i=0;i<n;i++) 
                nans[i]=-1; //nans用于记录列标，预置为-1 
        for(k=0;k<n;k++) //按列标循环 
        { 
                bigst=-1; 
                for(i=0;i<n;i++) //按行标循环 
                { 
                        if(nans[i]<0) 
                        //如果nans[i]>0说明第i行已作为某列的主元,不能再做其它列的主元 
                        { 
                                if(bigst-fabs(dMatrix[i][k])<0) //寻找某一列中的最大值作主元 
                                { 
                                        m=i; //m记录最大值所在的列数 
                                        bigst=fabs(dMatrix[i][k]); 
                                } 
                        } 
                } 
                if(bigst-min<=0) 
                //如果某列所有值都小于min，高斯消除法将引起较大误差,求解无意义 
                { 
                        nstop=-1; 
                        return(nstop); 
                } 
                pivot=dMatrix[m][k]; //用最大值作为消除某一列的主元，这将减小误差 
                nans[m]=k; //第k列的主元在第m行,或者说m行是k列的主元 
                w=-1/pivot; //消除因子 
                if(k<n-1) 
                { 
                        for(j=k+1;j<n;j++) //消去k列 
                        { 
                                v=dMatrix[m][j]*w; 
                                if(fabs(v)-epss>0) 
                                { 
                                        for(i=0;i<n;i++) //将m行的倍数加到其它行 
                                                dMatrix[i][j]=dMatrix[i][j]+dMatrix[i][k]*v; 
                                } 
                                dMatrix[m][j]=-v; 
                        } 
                } 
                v=an[m]*w; 
                if(fabs(v)-epss>0) 
                { 
                        for(i=0;i<n;i++) //增广矩阵常数项随以上的行变换而变 
                                an[i]=an[i]+dMatrix[i][k]*v; 
                } 
                an[m]=-v; 
        } 
        for(i=0;i<n;i++) 
                dMatrix[i][0]=an[i]; 
        for(i=0;i<n;i++) 
                an[nans[i]]=dMatrix[i][0]; 
        return(nstop); 
}