高斯算法顯示意外的行爲

Question

我已經實現了無樞軸的高斯算法。

import matplotlib.pyplot as plt 
import numpy as np 
import scipy as sp 

def gauss_solve(A,b): 
    """ 
    args: coefficient matrix A of dim(nxn) and vector b of dim(n) 
    of a system of linear equations with n unknowns. 
    note: no zeroes on the main diagonal of A allowed! 

    returns: vector x of dim(n) which solves the SLE 
    """ 
    while np.ndim(A) != 2 or A.shape[0] != A.shape[1]: 
     A = input(["The matrix you entered is not square, specify new input matrix A: "]) 
# print "A ok." 
    while np.ndim(b) != 1 or A.shape[1] != b.shape[0]: 
     b = input(["The dimension of the constant vector b is incorrect, please specify new input vector b"]) 
# print "b ok." 
    if np.linalg.det(A) == 0: 
     return "This linear system doesn't have a single unique solution." 
# print "System does have solution: " 
    n = len(b) 
    for i in xrange(n): # create triangular matrix 
     if A[i,i] == 0: 
      return "This implementation doesn't allow A to have zero entries on the main diagonal." 
     A[i] = A[i]/float(A[i,i]) 
     b[i] = b[i]/float(A[i,i]) 
     for l in xrange(i+1,n): 
      A[l] -= A[i]*A[l,i] 
      b[l] -= b[i]*A[l,i] 
    r = np.zeros(n) # result 
    for i in xrange(n): 
     r[-(i+1)] = b[-(i+1)] - np.dot(r,A[-(i+1)]) 
    return r 

def test_gauss(): 
    m = 10 
    e = 0.1 
    A = sp.rand(m,m) 
# A,b = np.array([[e,1.],[1.,1.]]),np.array([1.,e]) 
    b = sp.rand(m) 
    print gauss_solve(A,b) 
    print "Build-in function says: \n", np.linalg.solve(A,b) 

test_gauss() 

測試函數可以生成用於Ab和隨機條目。一切都運行完美我想，但我這裏有一個矩陣導致意想不到的結果：

A = [[e 1] [1 1]] 
b = [1 e] 

對於e != 1的解析解

x = [-1 e+1] 

但我試過e一些值，我只是不」得到分析解決方案。即使內置函數solve(A,b)失敗。例如x的第一個條目總是0（儘管它應該是-1，完全獨立於e）。有誰能解釋爲什麼會發生這種情況？

Answer 1

您的並行更新A和b是不正確的，因爲你正在使用新值A更新b。您需要更換線路：

A[i] = A[i]/float(A[i,i]) 
b[i] = b[i]/float(A[i,i]) 

的東西，如：

divisor = A[i,i] 
A[i] = A[i]/float(divisor) 
b[i] = b[i]/float(divisor) 

同樣，該行：

​​

與

multiplier = A[l,i] 
A[l] -= A[i]*multiplier 
b[l] -= b[i]*multiplier 

在你原來的代碼，行號爲b什麼都不做（不考慮浮點精度問題）：代碼的第一部分將b[i]除以1.0，而第二部分從b[l]減去0.0時間b[i]。