解决与gausian emlinination与蟒蛇枢轴的9x9矩阵

Question

我需要解决整个范围的8x8和9x9矩阵，所以我认为我可以建立一个Python程序，使整个事情更容易。

到目前为止我已成功地创建：

from __future__ import division 
import numpy as np 

def solveEqns(A,v): 
    def lu(A): 
     #Factor A into LU by Gaussian elimination with scaled partial pivoting 
     n, m = np.shape(A) 
     if n != m: 
      print "Error: input matrix is not square" 
      return None 
      # Generate initial index vector 
     p = range(n) 
      # Determine the largest (in magnitude) element in each row. These 
      # factors are used to scale the pivot elements for comparison purposes 
      # when deciding which row to use as a pivot row. 
     s = [0] * n 
     for i in xrange(n): 
      smax = 0.0 
      for j in xrange(n): 
       smax = max(smax, abs(A[i][j])) 
      s[i] = smax 
      # Begin Gaussian elimination. 
     for k in xrange(n - 1): 
     # Find the remaining row with the largest scaled pivot. 
      rmax = 0.0 
      for i in xrange(k, n): 
       r = abs(A[p[i][k]]/s[p[i]]) 
       if r > rmax: 
        rmax = r 
        j = i 
     # Row j has the largest scaled pivot, so "swap" that row with the 
     # current row (row k). The swap is not actually done by copying rows, 
     # but by swaping two entries in an index vector. 
      p[j], p[k] = (p[k], p[j]) 
     # Now carry out the next elimination step as usual, except for the 
     # added complication of the index vector. 
      for i in xrange(k + 1, n): 
       xmult = A[p[i],k]/A[p[k],k] 
       A[p[i],k] = xmult 
       for j in xrange(k + 1, n): 
        A[p[i],j] = A[p[i],j] - xmult * A[p[k],j] 
    # All done, return factored matrix A and permutation vector p 
     return (A, p) 

def solve(A, p, b): 
#Solves Ax = b given an LU factored matrix A and permuation vector p 
    n, m = np.shape(A) 
    if n != m: 
     print "Error: input matrix is not square" 
     return None 
# Forward solve 
    x = np.zeros(n) 
    for k in xrange(n - 1): 
     for i in xrange(k + 1, n): 
      b[p[i]] = b[p[i]] - A[p[i],k] * b[p[k]] 
# Backward solve 
    for i in xrange(n - 1, -1, -1): 
     sum = b[p[i]] 
     for j in xrange(i + 1, n): 
      sum = sum - A[p[i],j] * x[j] 
      x[i] = sum/A[p[i],i] 

# All done, return solution vector 
    return x 

lu(A) 
return solve(A,p,v) 

DEF电路（）：

A = np.array([[1,0,0,0,0,8,0,0,0],[0,1,0,0,5,0,0,0,0],[0,1,0,0,5,0,0,0,0],[0,0,0,1,-1,1,0,0,0],[0,0,1,0,0,0,1,-1,0],[0,0,1,0,0,0,1,0,-1],[0,1,0,0,-1,0,0,0,1],[1,0,0,0,0,-1,0,1,0],[1,-1,0,1,0,0,0,0,0]]) 

v = np.array([9,-12,-0.5,0,0,0,0,0,0]) 
I = solveEqns(A,v) 
return I 

解决9x9的矩阵A在末端。这是我需要解决的更容易的一个，所以可以在python之外解决它，以检查结果是否准确。

即时得到上线26回溯错误：

回溯（最近通话最后一个）：

File "<ipython-input-110-6daf773db1e3>", line 1, in <module> 
    solveEqns(A,b) 

    File "C:/Users/SamMc/Documents/Python Scripts/q6u1510416 v4.py", line 65, in solveEqns 
    lu(A) 

    File "C:/Users/SamMc/Documents/Python Scripts/q6u1510416 v4.py", line 26, in lu 
    r = abs(A[p[i][k]]/s[p[i]]) 

TypeError: 'int' object has no attribute '__getitem__' 

我无法弄清楚为什么它不通过从矩阵的一些拉动。

任何帮助将不胜感激。

感谢

山姆

Answer 1

您可以使用通过缩放旋转高斯消元法。代码如下所示。

import numpy as np 
def gauss_pivot(a,b,tol=1.0e-12): 
    """ 
    x = gaussPivot(a,b,tol=1.0e-12). 
    Solves [a]{x} = {b} by Gauss elimination with 
    scaled row pivoting 
    """ 
    a = np.copy(a) 
    b = np.copy(b) 
    n = len(b) 
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix 
    # Set up scale factors 
    s = np.zeros(n) 
    for i in range(n): 
     s[i] = max(np.abs(a[i,:])) # find the max of each row 
    for k in range(0, n-1): #pivot row 
     # Row interchange, if needed 
     p = np.argmax(np.abs(a[k:n,k])/s[k:n]) # find which row has max item for each col k, and scale by s 
     if abs(a[p,k]) < tol: 
      raise Exception("Matrix is singular") 
     if p != k: # swap rows if current row does not contain max item with the one contains max item within same col 
      a[[k,p+k],:] = a[[p+k, k],:] 
      b[k],b[p+k] = b[p+k],b[k] 
      s[k],s[p+k] = s[p+k],s[k] 

     # Elimination phase of matrix a 
     for i in range(k+1,n): 
      if a[i,k] != 0.0: # skip if a(i,k) is already zero 
       lam = a [i,k]/a[k,k] 
       a[i,k:n] = a[i,k:n] - lam*a[k,k:n] 
       b[i] = b[i] - lam*b[k] 

    if abs(a[n-1,n-1]) < tol: 
     raise Exception("Matrix is singular") 

    # Back substitution phase, solution is substituted by b 
    x = np.zeros_like(b) 
    x[n-1] = b[n-1]/a[n-1,n-1] 
    for k in range(n-2,-1,-1): 
     x[k] = (b[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k] 
    return x 

a = np.random.randn(100,100)*10 
b = np.random.randn(100)*10 

x = gauss_pivot(a,b) 

if np.allclose(np.dot(a,x), b) == True: 
    print("x is the correct solution") 

如果您希望代码执行得更快，您可能会用b替换x，所以在函数返回b时包含解决方案。 您可能还会略微修改消除阶段，因此矩阵a的一个对角线以下的元素不会归零，因为在后置替换阶段无关紧要。因此，代码变，如下所示：

import numpy as np 
def gauss_pivot(a,b,tol=1.0e-12): 
    """ 
    x = gaussPivot(a,b,tol=1.0e-12). 
    Solves [a]{x} = {b} by Gauss elimination with 
    scaled row pivoting 
    """ 
    a = np.copy(a) 
    b = np.copy(b) 
    n = len(b) 
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix 
    # Set up scale factors 
    s = np.zeros(n) 
    for i in range(n): 
     s[i] = max(np.abs(a[i,:])) # find the max of each row 
    for k in range(0, n-1): #pivot row 
     # Row interchange, if needed 
     p = np.argmax(np.abs(a[k:n,k])/s[k:n]) # find which row has max item for each col k, and scale by s 
     if abs(a[p,k]) < tol: 
      raise Exception("Matrix is singular") 
     if p != k: # swap rows if current row does not contain max item with the one contains max item within same col 
      a[[k,p+k],:] = a[[p+k, k],:] 
      b[k],b[p+k] = b[p+k],b[k] 
      s[k],s[p+k] = s[p+k],s[k] 

     # Elimination phase of matrix a 
     for i in range(k+1,n): 
      if a[i,k] != 0.0: # skip if a(i,k) is already zero 
       lam = a [i,k]/a[k,k] 
       a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n] 
       b[i] = b[i] - lam*b[k] 

    if abs(a[n-1,n-1]) < tol: 
     raise Exception("Matrix is singular") 

    # Back substitution phase, solution is substituted by b 

    b[n-1] = b[n-1]/a[n-1,n-1] 
    for k in range(n-2,-1,-1): 
     b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k] 
    return b 

要使用LU分解，而不是这是更理想的为含一个以上的列B中，LU代码如下所示

import numpy as np 
def lu_decomp(a,tol=1.0e-9): 
    a = np.copy(a) 
    n = len(a) 
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix 
    seq = np.arange(n, dtype=int) 
    s = np.zeros((n)) 
    for i in range(n): 
     s[i] = max(abs(a[i,:])) 
    for k in range(0,n-1): 
     p = np.argmax(np.abs(a[k:n,k])/s[k:n]) 
     if abs(a[p,k]) < tol: 
      raise Exception("Matrix is singular") 
     if p != k: 
      a[[k,p+k],:]  = a[[p+k, k],:] 
      s[k],s[p+k]  = s[p+k],s[k] 
      seq[k], seq[p+k] = seq[p+k],seq[k] 

     # Elimination 
     for i in range(k+1,n): 
      if a[i,k] != 0.0: 
       lam = a[i,k]/a[k,k] 
       a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n] 
       a[i,k] = lam 
    return a,seq 

def lu_solve(a,b,seq): 
    n = len(a) 
    x = b.copy() 
    for i in range(n): 
     x[i] = b[seq[i]] 
    # Solution 
    for k in range(1,n): 
     x[k] = x[k] - np.dot(a[k,0:k],x[0:k]) 
    x[n-1] = x[n-1]/a[n-1,n-1] 
    for k in range(n-2,-1,-1): 
     x[k] = (x[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k] 
    return x 

a2 = np.random.randn(500,500)*100 
b2 = np.random.randn(500,20)*100 

a_decomposed, seq = lu_decomp(a2) 

x2 = np.zeros_like(b2) 

for col in range(b2.shape[1]): 
    x2[:,col] = lu_solve(a_decomposed, b2[:, col], seq) 


if np.allclose(np.dot(a2,x2), b2) == True: 
    print("x2 is the correct solution") 

两种方法给出了输出，

高斯消

x是正确的溶液

LU方法

X2是正确的解决方案

我建议你使用SciPy的linalg包，从scipy.linalg进口解决，lu_factor，lu_solve。 他们执行方式更快的大矩阵大小。你可以使用上面相同的代码，但用numba jit注释它们，所以对于大型矩阵，性能会更好。

from numba import jit 
@jit 
def gauss_pivot(a, b): 
    ... 
    ... 

确认：由教授在科学和工程账面数值方法与Python代码的启发Jaan Kiusalaas

https://www.amazon.co.uk/Numerical-Methods-Engineering-Python-3/dp/1107033853/ref=sr_1_1?ie=UTF8&qid=1517845946&sr=8-1&keywords=numerical+method+in+science+and+engineering+with+python