二次项最小二乘法拟合代码

自网上一组代码的修改.

void Matrix2Order(double *dx, double *dy, int n, double* pCoef)

{

    int i = 0 , j = 0 , k = 0 ;

    double dMatrix[4][5];

    double * dTemp= new double [n];    //store sigma(xi*xi*xi)

    dMatrix[1][4]=0;

    for ( j = 0 ; j < n ; j ++ )

    {

        dTemp[j] = 1;

        dMatrix[1][4] += dy[j];

    }

    for ( i = 1 ; i < 3 ; i++)

    {

        dMatrix[1][i+1] = 0;

        dMatrix[i+1][4] = 0;

        for ( j = 0 ; j < n ; j++ )

        {

            dTemp[j] *= dx[j];

            dMatrix[1][i+1] += dTemp[j];

            dMatrix[i+1][4] += dTemp[j] * dy[j];

        }

    }

    for ( ; i < 5 ; i ++ )

    {

        dMatrix[3][i-1] = 0;

        for ( j = 0 ; j < n ; j++ )

        {

            dTemp[j] *= dx[j];

            dMatrix[3][i-1] += dTemp[j];

        }

    }

    dMatrix[1][1]=n;

    dMatrix[3][1]=dMatrix[1][3];

    dMatrix[2][1]=dMatrix[1][2];

    dMatrix[2][2]=dMatrix[1][3];

    dMatrix[2][3]=dMatrix[3][2];

    delete dTemp;

    //calculate the three variables equation

    for( k=1;k<3;k++) //消元过程

    {

        for( i=k+1;i<4;i++)

        {

            double p1=0;

            if(dMatrix[k][k]!=0)

                p1=dMatrix[i][k]/dMatrix[k][k];

            for(int j=k;j<5;j++)

                dMatrix[i][j]=dMatrix[i][j]-dMatrix[k][j]*p1;

        }

    }

    pCoef[3]=dMatrix[3][4]/dMatrix[3][3];

    for(int l=2;l>=1;l--)//回代求解

    {

        double sum=0;

        for(i=l+1;i<=3;i++)

            sum+=dMatrix[l][i]*coefficient[i];

        pCoef[l]=(dMatrix[l][4]-sum)/dMatrix[l][l];

    }

}