天天看點

高斯消元法和列主銷元法(Python實作)

目錄

​​1、原理​​

​​(1)思維導圖​​

​​ (2)原理​​

​​2、案例及實作 ​​

​​(1)案例​​

​​(2)代碼實作 ​​

​​(3)結果​​

​​3、誤差分析及心得體會​​

1、原理

(1)思維導圖

高斯消元法和列主銷元法(Python實作)

 (2)原理

高斯消元法和列主銷元法(Python實作)

2、案例及實作 

(1)案例

(2)代碼實作 

import numpy as np
# 判斷系數矩陣是否為非奇異矩陣
def IsItNonSingular(A):
    row = len(A)  # 系數矩陣行數
    col = len(A[0])  # 系數矩陣列數
    # 控制第i步,A矩陣為n階方陣
    for i in range(row - 1):
        # 若主元元素為0,則進行行變換
        if A[i][i] == 0:
            # 從後往前找出主元為0的列中,首個不為零的行
            # 對主元為0的行與此行進行行變換
            Str = []
            for h in range(row - 1, i, -1):  # 添加該列所有元素
                Str.append(A[h][i])
            Num = col - Str.index(next(filter(lambda x: x != 0, Str)))
            A[Num - 1], A[i] = A[i], A[Num - 1]  # 做行變換
        # 消元運算
        for j in range(i + 1, row, 1):
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1):
                A[j][k] = A[j][k] - coeff * A[i][k]  # 系數矩陣消元
    # 判斷該方程組系數矩陣是否為非奇異矩陣
    for i in range(row):
        if A[i][i] == 0:
            print("Coefficient matrix is not a nonsingular matrix.")
            return 'N'
    print("Coefficient matrix is a nonsingular matrix.")


def GuassElimination(A, b):
    row = len(A)  # 系數矩陣行數
    col = len(A[0])  # 系數矩陣列數
    ε = 1E-5  # 定義一個小量
    print("coefficient matrix:", A)
    print("Constant column:", b)
    # 控制第i步,高斯消元需要n-1步,A矩陣為n階方陣
    for i in range(row - 1):
        if abs(A[i][i]) <= ε:  # 若主元為一個小量,則采用列主消元
            return None
        # 消元運算
        for j in range(i + 1, row, 1):
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1):
                A[j][k] = A[j][k] - coeff * A[i][k]  # 系數矩陣消元
            b[j] = b[j] - coeff * b[i]  # 對應常數列消元

    # 回代過程
    x = [0] * col  # 初始化元組,用于後面存放解
    x[col - 1] = b[col - 1] / A[col - 1][col - 1]  # 第n個解
    for i in range(row - 2, -1, -1):
        for j in range(col - 1, i, -1):
            b[i] = b[i] - A[i][j] * x[j]
        x[i] = b[i] / A[i][i]
    print("The solution of the equations is:")
    for i in range(col):
        print("x{}:".format(i + 1), '%.8f' % x[i])
    print("\n")
    return x


def ColumnPrincipalElimination(A, b):
    row = len(A)  # 系數矩陣行數
    col = len(A[0])  # 系數矩陣列數
    print("coefficient matrix:", A)
    print("Constant column:", b)
    # 控制第i步,消元需要n-1步,A矩陣為n階方陣
    for i in range(row - 1):
        # 每步運算前找出列中絕對值最大元素
        # 作行變換,讓絕對值最大的元素行作主元
        Str = []
        for j in range(i, row, 1):
            Str.append(A[j][i])
        Num = Str.index(max(Str)) + i
        A[Num], A[i] = A[i], A[Num]  # 行變換
        b[Num], b[i] = b[i], b[Num]

        # 消元運算
        for j in range(i + 1, row, 1):
            coeff = A[j][i] / A[i][i]
            for k in range(i, col, 1):
                A[j][k] = A[j][k] - coeff * A[i][k]  # 系數矩陣消元
            b[j] = b[j] - coeff * b[i]  # 對應常數列消元

    # 回代過程
    x = [0] * col  # 初始化元組,用于後面存放解
    x[col - 1] = b[col - 1] / A[col - 1][col - 1]  # 第n個解
    for i in range(row - 2, -1, -1):
        for j in range(col - 1, i, -1):
            b[i] = b[i] - A[i][j] * x[j]
        x[i] = b[i] / A[i][i]
    print("The solution of the equations is:")
    for i in range(col):
        print("x{}:".format(i + 1), '%.8f' % x[i])
    # print("\n")
    return x
def main():
    b1 = [0.4043, 0.1550, 0.4240, -0.2557]
    b2 = [0.4043, 0.1550, 0.4240, -0.2557]
    test1 = [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.1920, 0.3781, 0.0643],
                  [0.1784, 0.4002, 0.2785, 0.3927]]
    test2 = [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.1920, 0.3781, 0.0643],
                  [0.1784, 0.4002, 0.2785, 0.3927]]
    # python為動态語言,
    # 在判斷非奇異矩陣過程中會改動初值
    # 故重新指派
    IsItNonSingular(test1)
    print("test1 Gaussian elimination test:")
    test1 = [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.1920, 0.3781, 0.0643],
                  [0.1784, 0.4002, 0.2785, 0.3927]]
    GuassElimination(test1, b1)

    IsItNonSingular(test2)
    print("test2Gaussian column principal elimination test:")
    test2 = [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.1920, 0.3781, 0.0643],
                  [0.1784, 0.4002, 0.2785, 0.3927]]
    ColumnPrincipalElimination(test2, b2)


if __name__ == '__main__':
    main()      

(3)結果

Coefficient matrix is a nonsingular matrix.
test1 Gaussian elimination test:
coefficient matrix: [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.192, 0.3781, 0.0643], [0.1784, 0.4002, 0.2785, 0.3927]]
Constant column: [0.4043, 0.155, 0.424, -0.2557]
The solution of the equations is:
x1: -0.18034012
x2: -1.66163443
x3: 2.21499710
x4: -0.44669701


Coefficient matrix is a nonsingular matrix.
test2Gaussian column principal elimination test:
coefficient matrix: [[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.192, 0.3781, 0.0643], [0.1784, 0.4002, 0.2785, 0.3927]]
Constant column: [0.4043, 0.155, 0.424, -0.2557]
The solution of the equations is:
x1: -0.18034012
x2: -1.66163443
x3: 2.21499710
x4: -0.44669701

Process finished with exit code 0      

3、誤差分析及心得體會