Gauss_Jordan消元法之Python实现
最近在看MIT OCW的线性代数,本来想自己写个Gauss_Jordan消元法的程序,结果前几天上维基百科的时候碰巧看到了,拷贝过来(Python2): [python] def gauss_jordan(m, eps = 1.0/(10**10)): “”“Puts given matrix (2D array) into the Reduced Row Echelon Form. Returns True if successful, False if ‘m’ is singular. NOTE: make sure all the matrix items support fractions! Int matrix will NOT work! Written by J. Elonen in April 2005, released into Public Domain”“” (h, w) = (len(m), len(m[0])) for y in range(0,h): maxrow = y for y2 in range(y+1, h): # Find max pivot if abs(m[y2][y]) > abs(m[maxrow][y]): maxrow = y2 (m[y], m[maxrow]) = (m[maxrow], m[y]) if abs(m[y][y]) <= eps: # Singular? return False for y2 in range(y+1, h): # Eliminate column y c = m[y2][y] / m[y][y] for x in range(y, w): m[y2][x] -= m[y][x] * c for y in range(h-1, 0-1, -1): # Backsubstitute c = m[y][y] for y2 in range(0,y): for x in range(w-1, y-1, -1): m[y2][x] -= m[y][x] * m[y2][y] / c m[y][y] /= c for x in range(h, w): # Normalize row y m[y][x] /= c return True [/python] 需要注意的是,由于Python2中整数除法返回整数,所以对于整数矩阵该程序不一定能正常工作,解决办法有:使用浮点数,或者改用Python3。因为Python3中/操作符固定返回浮点数。 下面有个调用的示例: [python] mat=[[1.0,2.0,1.0,2.0],[3.0,8.0,1.0,12.0],[0.0,4.0,1.0,2.0]] print mat gauss_jordan(mat) print mat [/python]