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

自网上一组代码的修改.

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];
    }
}