例題6-1 利用高斯消去法求線性代數的值

'''
例題6-1
利用高斯消去法求線性代數的值
E1 :  -1 x1 +    x2  + 2 x3 = 2
E2 :   3 x1  - 1 x2  +    x3 = 6

E3 :  -1 x1 + 3 x2  + 4 x3 = 4

/* Based on Gaussian Elimination method
 * for solving the n x n linear algebra system
 * a11 x1+a12 x2+...+a1n xn=b1
 * a21 x1+a22 x2+...+a2n xn=b2
 * .       .         .       .
 * .       .         .       .
 * an1 x1+an2 x2+...+ann xn=bn
 * Input number of unknowns and equations n
 * with coefficent a11,a12,...,ann and b1,b2,
 * ...bn. Output solution x1,x2,x3,...,xn.
 */
 '''

def pprint(A):
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n+1):
            line += str(A[i][j]) + "\t"
            if j == n-1:
                line += "| "
        print(line)
    print("")


def gauss(A):
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n+1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i+1, n):
            c = -A[k][i]/A[i][i]
            for j in range(i, n+1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = A[i][n]/A[i][i]
        for k in range(i-1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x
#============================================
A=[[-1.0 ,  1.0 , 2.0 , 2.0],
      [3.0 , -1.0 , 1.0,  6.0],
      [-1.0 ,  3.0 , 4.0 , 4.0]]
n=3
# Print input
pprint(A)
# Calculate solution
x = gauss(A)

# Print result
print ("\nThe values of the variables are as follows:\n")
for i in range(0 ,n):
        print( "x",i+1,"=",round(x[i],4))         #// Print the values of x, y,z,....   

輸出畫面

======== RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX6-1-3.py ===========

-1.0 1.0 2.0 | 2.0

3.0 -1.0 1.0 | 6.0

-1.0 3.0 4.0 | 4.0

The values of the variables are as follows:

x 1 = 1.0

x 2 = -1.0

x 3 = 2.0

>>>