matlab约当消去法,Gauss-Jordan 高斯约当消去法

我听取了我们老师的想法，把原书中的命名改为以意义为中心的，不简化的命名方式。原本函数名为gaussj，现在改为gaussJordan

Gauss-Jordan消去法解的是这个线性方程组集合：(我完全不会LaTex呢，而且鉴于复制贴贴的方便，或者转载分享之后会出现不能显示的问题，我会尽量少用Latex公式的)

[tex]A*(X_0 \cup X_1 \cup X_2 \cup Y) = (b_0 \cup b_1 \cup b_2 \cup 1)[/tex]

如果上面一行的公式显示不出，请看这行：A * (X0 || X1 || X2 || Y) = (b0 || b1 || b2 || 1)

其中A是一个n*n的非奇异矩阵，b是n*m的矩阵。

整个式子相当于拆分成4个子式子：

A * X0 = b0

A * X1 = b1

A * X2 = b2

A * Y = 1

算法的运算思想是：

施行初等变换将A变为单位矩阵，则

Xi = A的逆矩阵 * b

Y = A的逆矩阵

算法的步骤：

(1)对于每一列k：

(2)选主元pivot element(绝对值最大的那个)(为了减小舍入误差)，记为a[i][j]

(3)如果选出的主元不在对角线上，我们把它换到对角线上

(4)如果主元在对角线上，做消元。

(5)最后再把交换过的元素交换回来以得到正确的A的逆矩阵。

/***@file gaussJordan.h*/#ifndef GAUSSJORDAN_H_INCLUDED#define GAUSSJORDAN_H_INCLUDED#include "nr3.h"/***@brief Solve linear equations A*X = b using Gauss-Jordan elimination.*@parama[0..n-1][0..n-1]Input as Matrix A in A * x = b, Output its matrix inverse*@paramb[0..n-1][0..m-1]Input as Matrix b in A * x = b, Output the corresponding set of solution vectors*Linear equation solution by Gauss-Jordan elimination, equation example as follow,*A * (X0 || X1 || X2 || Y) = (b0 || b1 || b2 || 1):**(a00 a01 a02 a03)((x00)(x01)(x02)(y00 y01 y02 y03))*(a10 a11 a12 a13)((x10)(x11)(x12)(y10 y11 y12 y13))*(a20 a21 a22 a23)*((x20)||(x21)||(x22)||(y20 y21 y22 y23))*(a30 a31 a32 a33)((x30)(x31)(x32)(y30 y31 y32 y33))**((b00)(b01)(b02)(1 0 0 1))*((b10)(b11)(b12)(0 1 0 0))*=((b20)||(b21)||(b22)||(0 0 1 0))*((b30)(b31)(b32)(0 0 0 1))**@attention A and Y should be square matrixes; b and X should be the column-vector matrix.*This function simultaneously solves the linear sets:*A * x0 = b0*A * x1 = b1*A * x2 = b2*A * Y = 1**Input matrix:*a[0..n-1][0..n-1], b[0..n-1][0..m-1] is input containing*the m right-hand side vectors*Output matrix:*a is replaced by its matrix inverse*b is replaced by the corresponding set of solution vectors*/void gaussJordan(MatDoub_IO &a, MatDoub_IO &b){Int i, icol, irow, j, k ,l , ll, n = a.nrows(), m = b.ncols();Doub big, dum, pivinv;VecInt indxc(n), indxr(n), ipiv(n, 0); //These integer arrays are used for bookkeeping on the pivoting.for(i = 0; i < n; ++i){//Main loop over the columns to be reducedbig = 0.0;for(j = 0; j < n; ++j){//Outer loop of the search for a pivot element.if(ipiv[j] != 1){for(k = 0; k < n; ++k){if(ipiv[k] == 0){if(abs(a[j][k]) >= big){big = abs(a[j][k]);irow = j;icol = k;}//if}//if}//for}//if} //for++(ipiv[icol]);/**We now have the pivot element, so we interchange rows, if needed, to put the pivot element on the diagonal.*The columns are not physically interchanged, only reabeled: indxc[i], the column of the (i+1)-th pivot element,*is the (i+1)-th column that is reduced, while indxr[i] is the row in which that pivot element was originally located.*If indxr[i] != indxc[i], there is an implied column interchange. With this form of bookkeeping, the solution b's will*end up in the correct order, and the inverse matrix will be scrambled by columns.*/if(irow != icol){for(l = 0; l < n; ++l)SWAP(a[irow][l], a[icol][l]);for(l = 0; l < m; ++l)SWAP(b[irow][l], b[icol][l]);} //ifindxr[i] = irow; //We are now readey to divide the pivot row by the pivot element, located at irow and icol.indxc[i] = icol;if(a[icol][icol] == 0.0)NRthrow("gaussj: Singular Matrix");pivinv = 1.0/a[icol][icol];a[icol][icol] = 1.0;for(l = 0; l < n; ++l)a[icol][l] *= pivinv;for(l = 0; l < m; ++l)b[icol][l] *= pivinv;for(ll = 0; ll < n; ++ll){//Next, we reduce the rowsif(ll != icol){//except for the pivot one, of course.dum = a[ll][icol];a[ll][icol] = 0.0;for(l = 0; l < n; ++l)a[ll][l] -= a[icol][l] * dum;for(l = 0; l < m; ++l)b[ll][l] -= b[icol][l] * dum;} //if} //for} //for/**This is the end of the main loop over columns of the reduction. It only remains to unscramble*the solution in view of the column interchanges. We do this by interchanging pairs of columns*in the reverse order that the permutation was built up.*/for(l = n-1; l >= 0; --l){if(indxr[l] != indxc[l])for(k = 0; k < n; ++k)SWAP(a[k][indxr[l]], a[k][indxc[l]]);} //for... And we are done.} //gaussj/***@brief Overloaded version with no right-hand sides. Replaces a by its inverse.*@parama[0..n-1][0..n-1]Input as Matrix A in A * x = 0, Output its matrix inverse*/void gaussJordan(MatDoub_IO &a){MatDoub_IO b(a.nrows(), 1, 0.0); //Dummy vector with zero columns.gaussJordan(a,b);}#endif // GAUSSJ_H_INCLUDED

大家如果想要做测试的话，可以试试以下代码：

/***@file gaussJordan_test.cpp*/#include #include "nr3.h"#include "gaussJordan.h"using namespace std;int main(){MatDoub_IO a,b;int rows, cols;//INPUT DATA:cout << "Please input rows, cols and the matrix data for A: " << endl;cin >> rows >> cols;a.resize(rows, cols);readMatrix_console(a);cout << "Please input rows, cols and the matrix data for B: " << endl;cin >> rows >> cols;b.resize(rows, cols);readMatrix_console(b);//Do the Linear SolvegaussJordan(a, b);//PRINTcout << "Gives the Inverse of A: " << endl;printMatrix_console(a);cout << "Gives the solution matrix: " << endl;printMatrix_console(b);return 0;}

欢迎加入我爱机器学习QQ14群：336582044

微信扫一扫，关注我爱机器学习公众号