高斯列主元消去法

Gauss.py

#-*-coding: utf-8 -*-


EPS = 1.0e-4


def guass(A, b):
    """ 高斯消元法
         输入要求：
            A = [[x11, x12, x13], [x21, x22, x23], [x31, x32, x33]]
            b = [b1, b2, b3]
    """

    if not len(A) == len(b):
        print "Error: Illegal input."
        return None

    n = len(A)

    i = 0;
    for line in A:
        main_element_index = i                  # 计算主元
        main_element = line[main_element_index]

        for j in range(i+1, n):
            if abs(main_element) < abs(line[j]):
                main_element = line[j]
                main_element_index = j
        
        if (abs(A[main_element_index][i]) < EPS):
            print "方程组有无数解！"
            return None

        if not i == main_element_index:
            swap_A(A, i, main_element_index)    # 换行
            swap_b(b, i, main_element_index)

        for j in range(i + 1, n):               # 消元
            m = A[j][i] / A[i][i]
            for k in range(i, n):
                A[j][k] = A[j][k] - m * A[i][k]
            b[j] = b[j] - m * b[i]

        i = i + 1

    if A[n - 1][n - 1] == 0:
        print "矩阵行列式值为0"
        return None

    x = b

    x[n - 1] = x[n - 1] / A[n - 1][n - 1]

    for i in range(n - 2 , -1, -1):
        for j in range(n - 1, i, -1):
            x[i] = x[i] - A[i][j] * x[j]
        x[i] = x[i] / A[i][i]

    return x


def swap_A(A, one, two):
    if not len(A[one]) == len(A[two]):
        print "Cannot swap two lines. Not equal length of them. "
        return None
    for i in range(0, len(A[one])):
        temp = A[two][i]
        A[two][i] = A[one][i]
        A[one][i] = temp


def swap_b(b, one, two):
    temp = b[one]
    b[one] = b[two]
    b[two] = temp


"""Example"""
A = [[1, -0.4, -0.3],
     [-0.3, 0.9, -0.2],
     [-0.2, -0.1, 0.8]]

b = [60000, 30000, 50000]

x = guass(A, b)

print x