求n阶矩阵的逆

 

多情况下我们需要用到矩阵，并且需要使用它的逆，下面介绍一下一种求矩阵逆的方法及实现，该方法取自《运筹学（第三版）》，有兴趣大家可以参考一下，若有错误，也希望大家指正。

设矩阵

A=

a11	a12	…	a1m
a21	a22	…	a2m
…	…	…	…
am1	am2	…	amm

求其逆矩阵可以从第一列开始，

P₁=

a11

a21

…

am1

以a11为主元变换为

1/a11

-a21/a11

…

-am1/a11

,然后构造含有该列但其他列均为单位列的矩阵

E₁=

1/a11	0	…	0
-a21/a11	1	…	0
…	…	…	…
-am1/a11	0	…	1

则有E₁P₁=

1

0

…

0

且E₁A=

1	a(1)12	…	a(1)1m
0	a(1)22	…	a(1)2m
…	…	…	…
0	a(1)m2	…	a(1)mm

然后再以第二列的a⁽¹⁾₂₂ 为主元变换为

-a(1)12/a(1)22

1/a(1)22

…

-a(1)m2/a(1)22

然后如同第一步一样构造E₂，计算出E₂E₁A，此时第一二列成为单位列，如此一步步进行，直到E_m…E₂E₁A为单位矩阵

可求得逆矩阵A^-1为E_m…E₂E₁  

 下面即是我自己实现的方法，使用了面向对象的方法：

 

 

#ifndef MaxArrayVar #define MaxArrayVar 30 #endif class Matrix { public: virtual int Run(); Matrix(); Matrix(float (*Array)[MaxArrayVar],unsigned int row); int ReSet(float (*Array)[MaxArrayVar],unsigned int row); virtual ~Matrix(); protected: virtual void OriginalMatrix(bool printOrNot); virtual void PrintMatrix(bool printOrNot=true); bool Enable; //矩阵是否可运算(是否初始化或包含可计算数据) virtual int Compute_Ei(float *Pi,unsigned int row,float *Ei); void Initialize(); virtual int Multiplicate(float (*Array1)[MaxArrayVar],float (*Array2)[MaxArrayVar],unsigned int row1,unsigned int common,unsigned int column2,float (*result)[MaxArrayVar]); virtual int Inverse(unsigned int column); unsigned int _row; //矩阵的行列数 //unsigned int _column; float _Matrix[MaxArrayVar][MaxArrayVar]; //初始矩阵 float _MatrixE[MaxArrayVar][MaxArrayVar]; //Em...E1 float _MatrixEA[MaxArrayVar][MaxArrayVar]; //Em...E1A bool printOrNot; //是否输出到屏幕 };

/*初始化矩阵,并获得E1,E1A*/ void Matrix::Initialize() { float P1[MaxArrayVar]={0},E1[MaxArrayVar]={0}; unsigned int rowIndex=0; if(Enable) { for(;rowIndex<_row;rowIndex++) { P1[rowIndex]=_Matrix[rowIndex][0]; } Compute_Ei(P1,0,E1); for(rowIndex=0;rowIndex<_row;rowIndex++) { _MatrixE[rowIndex][0]=E1[rowIndex]; } Multiplicate(_MatrixE,_Matrix,_row,_row,_row,_MatrixEA); } } /*重置矩阵,并初始化*/ int Matrix::ReSet(float (*Array)[MaxArrayVar],unsigned int row) { unsigned int rowIndex=0; if(row>MaxArrayVar) { Enable=false; _row=0; return -1; } for(;rowIndex<row;rowIndex++) { _MatrixE[rowIndex][rowIndex]=1; for(unsigned int column=0;column<row;column++) { _Matrix[rowIndex][column]=Array[rowIndex][column]; _MatrixEA[rowIndex][column]=0; } } _row=row; Enable=true; Initialize(); return 0; } /*矩阵乘法运算*/ int Matrix::Multiplicate(float (*Array1)[MaxArrayVar],float (*Array2)[MaxArrayVar],unsigned int row1,unsigned int common,unsigned int column2,float (*result)[MaxArrayVar]) {/*Debug*/ float tmpresult[MaxArrayVar][MaxArrayVar]={0}; if(MaxArrayVar<=row1||MaxArrayVar<=common||MaxArrayVar<=column2) return -1; unsigned int rowIndex=0,columnIndex=0; …… return 0; }/*Debug*/ /*求解ξi*/ int Matrix::Compute_Ei(float *Pi, unsigned int row,float *Ei) { unsigned int rowIndex=0; float Ai=Pi[row]; if(Ai==0) return -1; Ai=-1/Ai; for(;rowIndex<_row;rowIndex++) { if(rowIndex!=row) { Ei[rowIndex]=Ai*Pi[rowIndex]; } else { Ei[rowIndex]=-Ai; } } return 0; } /*求解Em…E1,Em…E1A,即迭代处理过程*/ int Matrix::Inverse(unsigned int column) {/*Debug*/ float Pi[MaxArrayVar]={0},Ei[MaxArrayVar]={0}; float MatrixEi[MaxArrayVar][MaxArrayVar]={0}; //float MatrixEA[MaxArrayVar][MaxArrayVar]={0}; unsigned int rowIndex=0; if(Enable) { for(;rowIndex<_row;rowIndex++) { Pi[rowIndex]=_MatrixEA[rowIndex][column]; MatrixEi[rowIndex][rowIndex]=1; } if(-1==Compute_Ei(Pi,column,Ei)) return -1; for(rowIndex=0;rowIndex<_row;rowIndex++) { MatrixEi[rowIndex][column]=Ei[rowIndex]; } /*Debug*/ if(-1==Multiplicate(MatrixEi,_MatrixEA,_row,_row,_row,_MatrixEA)) return -1; if(-1==Multiplicate(MatrixEi,_MatrixE,_row,_row,_row,_MatrixE)) return -1; } return 0; } /*迭代运算,求解矩阵的逆*/ int Matrix::Run() {/*Debug*/ if(_row<1) return -1; if(_row<2) return 1; OriginalMatrix(printOrNot); PrintMatrix(printOrNot); for(unsigned int rowIndex=1;rowIndex<_row;rowIndex++) {/*Debug*/ if(-1==Inverse(rowIndex)) return -1; PrintMatrix(printOrNot); } return 0; }

学习运筹学时候实现的，在此贴出来，希望大家能与我交流（别扔砖头啊）^_^！