socketref,再见!高德

https://github.com/adoggie

  C++博客 :: 首页 :: 联系 :: 聚合  :: 管理
  246 Posts :: 4 Stories :: 312 Comments :: 0 Trackbacks

常用链接

留言簿(54)

我参与的团队

搜索

  •  

最新评论

阅读排行榜

评论排行榜

It can be used to solve linear equation systems or to invert a matrix.
高斯消元法用于解决线性代数求多元方程组的解,或者用于求可逆矩阵的逆
嘿嘿,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
使用方法 :

If your matrix is of form [A:x] (as is usual when solving systems), items of A and x both have to be divisible by items of A but not the other way around. Thus, you could, for example, use floats for A and vectors for x. Example:
mtx = [[1.0, 1.0, 1.0, Vec3(0.0,  4.0, 2.0), 2.0],
       [2.0, 1.0, 1.0, Vec3(1.0,  7.0, 3.0), 3.0],
       [1.0, 2.0, 1.0, Vec3(15.0, 2.0, 4.0), 4.0]]
if gauss_jordan(mtx):
  print mtx
else:
  print "Singular!"
# Prints out (approximately):
#
# [[1.0, 0.0, 0.0, (  1.0,  3.0,  1.0),  1.0],
#  [0.0, 1.0, 0.0, ( 15.0, -2.0,  2.0),  2.0],
#  [0.0, 0.0, 1.0, (-16.0,  3.0, -1.0), -1.0]]
Auxiliary functions contributed by Eric Atienza (also released in Public Domain):
def solve(M, b):
  """
  solves M*x = b
  return vector x so that M*x = b
  :param M: a matrix in the form of a list of list
  :param b: a vector in the form of a simple list of scalars
  """
  m2 = [row[:]+[right] for row,right in zip(M,b) ]
  return [row[-1] for row in m2] if gauss_jordan(m2) else None
def inv(M):
  """
  return the inv of the matrix M
  """
  #clone the matrix and append the identity matrix
  # [int(i==j) for j in range_M] is nothing but the i(th row of the identity matrix
  m2 = [row[:]+[int(i==j) for j in range(len(M) )] for i,row in enumerate(M) ]
  # extract the appended matrix (kind of m2[m:,...]
  return [row[len(M[0]):] for row in m2] if gauss_jordan(m2) else None
def zeros( s , zero=0):
    """
    return a matrix of size `size`
    :param size: a tuple containing dimensions of the matrix
    :param zero: the value to use to fill the matrix (by default it's zero )
    """
    return [zeros(s[1:] ) for i in range(s[0] ) ] if not len(s) else zero



算法伪代码:
i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot 
in column j, starting in row i:
  maxi :
= i
  
for k := i+1 to m do
    
if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi :
= k
    end 
if
  end 
for
  
if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but 
do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry 
in row i by A[i,j]
    Now A[i,j] will have the value 
1.
    
for u := i+1 to m do
      subtract A[u,j] 
* row i from row u
      Now A[u,j] will be 
0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end 
for
    i :
= i + 1
  end 
if
  j :
= j + 1
end 
while


posted on 2011-08-29 02:47 放屁阿狗 阅读(503) 评论(0)  编辑 收藏 引用 所属分类: Math

只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理