用高斯消元法求M行N列的方程(含N个未知数,N个方程)的解

来源：互联网发布：苹果手机铃声软件推荐编辑：程序博客网时间：2024/04/27 17:05

// 函数名: void GetCoeffEx(double **matrixTmp,const int iRow,const int iCol)

// 出口参数: 如果定义了一个 matrix[m][n] 则 matrix[1..m-1][n-1]的值就是相应的 a,b,c....的参数

// 下面的文件可以直接在C++ (VC6) 中执行。

// 作者: tonyjqq1 Email: tonyjqq@163.com (如转贴请加入此行)

#include <malloc.h>
#include <stdio.h>

void GetCoeffEx(double ** matrixTmp,const int iRow,const int iCol)
{

int i,j,k,z;

double** matrix;
matrix=(double **)malloc(sizeof(double *)*iRow);
for (i=0;i<iRow;i++)
matrix[i]=(double *)malloc(sizeof(double)*iCol);

double** matrix1;
matrix1=(double **)malloc(sizeof(double *)*iRow);
for (i=0;i<iRow;i++)
matrix1[i]=(double *)malloc(sizeof(double)*iCol);

for(i=0;i<iRow;i++)
for(j=0;j<iCol;j++)
     {
       matrix[i][j]    = *( (double*) matrixTmp +(i* iCol +j) );
      matrix1[i][j] = matrix[i][j] ;
     }

for(z=0;z<iRow;z++)
     {
       for(i=z;i<iRow;i++)
       for(j=z;j<iCol;j++)
    {
      for(k=z;k<iRow;k++)
         {
    if(k!=i) matrix[i][j] = matrix[i][j] * matrix1[k][z] ;
         }
   if(i!=z) matrix[i][j] = matrix[i][j] - matrix[z][j];
    }


    for(i=0;i<iRow;i++)
    for(j=0;j<iCol;j++)
    {
          matrix1[i][j] = matrix[i][j] ;
    }
     }

for (i=0;i<iRow;i++)
free(matrix1[i]);

free(matrix1);
matrix1 = NULL;

for(i=0;i<iRow;i++)
     {
    matrix[iRow-1-i][iCol-1] = matrix[iRow-1-i][iCol-1] / matrix[iRow-1-i][iCol-2-i];
       for(z = 0; z<(iRow-1-i); z++) matrix[z][iCol-1] = matrix[z][iCol-1] - matrix[z][iCol-2-i] * matrix[iRow-1-i][iCol-1];
       *((double*) matrixTmp + (iRow-i-1) * iCol + iCol-1 ) = matrix[iRow-1-i][iCol-1];

}

for (i=0;i<iRow;i++)
free(matrix[i]);

free(matrix);
matrix = NULL;

}

void GetCoeff(double matrix[3][4])
{
double matrix1[3][4];

int i,j;

for( i=0;i<=2;i++)
for( j=0;j<=3;j++)
matrix1[i][j] = matrix[i][j];

for( i=0;i<=3;i++)
{
  matrix[0][i] = matrix1[0][i] * matrix1[1][0] * matrix1[2][0];
  matrix[1][i] = matrix1[1][i] * matrix1[0][0] * matrix1[2][0];
  matrix[2][i] = matrix1[2][i] * matrix1[0][0] * matrix1[1][0];
}

for(i=0;i<=3;i++)
{
   matrix[1][i] = matrix[1][i] - matrix[0][i];
  matrix[2][i] = matrix[2][i] - matrix[0][i];
}

for( i=0;i<=2;i++)
for( j=0;j<=3;j++)
matrix1[i][j] = matrix[i][j];

for(j=1;j<=3;j++)
{
matrix[1][j] = matrix[1][j] * matrix1[2][1];
matrix[2][j] = matrix[2][j] * matrix1[1][1];
}

for(j=1;j<=3;j++) matrix[2][j] = matrix[2][j] - matrix[1][j];

matrix[2][3] = matrix[2][3]/ matrix[2][2];

matrix[1][3] = matrix[1][3] - matrix[2][3] * matrix[1][2];
matrix[0][3] = matrix[0][3] - matrix[2][3] * matrix[0][2];

matrix[1][3] = matrix[1][3] / matrix[1][1];

matrix[0][3] = matrix[0][3] - matrix[0][1] * matrix[1][3];

matrix[0][3] = matrix[0][3] / matrix[0][0];
}

void main()
{
double matrix[3][4];

    double NLCal = 489753.0, HLCal = 1373902.0 , FLCal = 2036800.0;
   double dCapacity= 35.000000;
    matrix[0][0]= NLCal * NLCal;
    matrix[0][1]= NLCal ;
    matrix[0][2]= 1.0;
    matrix[0][3]= 0.0;

    matrix[1][0]= HLCal * HLCal;
    matrix[1][1]= HLCal ;
    matrix[1][2]= 1.0;
    matrix[1][3]= 0.5 * dCapacity * 1000.0;

    matrix[2][0]= FLCal * FLCal;
    matrix[2][1]= FLCal ;
    matrix[2][2]= 1.0;
    matrix[2][3]= dCapacity * 1000.0;

  /*
matrix[0][0]= 934368.0 * 934368.0;
matrix[0][1]= 934368.0 ;
matrix[0][2]= 1.0;
matrix[0][3]= 0.0;

matrix[1][0]= 1589152.0 * 1589152.0;
matrix[1][1]= 1589152.0 ;
matrix[1][2]= 1.0;
matrix[1][3]= 3.0;

matrix[2][0]= 2243903.0 * 2243903.0;
matrix[2][1]= 2243903.0 ;
matrix[2][2]= 1.0;
matrix[2][3]= 6.0;
*/

//GetCoeff(matrix);

/*

matrix[0][0]= 1;
matrix[0][1]= 2;
matrix[0][2]= -3;
matrix[0][3]= 0;

matrix[1][0]= 2;
matrix[1][1]= 2;
matrix[1][2]= 4;
matrix[1][3]= 8;

matrix[2][0]= 3;
matrix[2][1]= 1 ;
matrix[2][2]= 2;
matrix[2][3]= 6;
*/

GetCoeffEx(reinterpret_cast<double **>(matrix),3,4);

printf("the value of a is: %e/n",matrix[0][3]);
printf("the value of b is: %e/n",matrix[1][3]);
printf("the value of c is: %e/n",matrix[2][3]);
}