python输出高斯消元法求解行列式的过程

前言

有人可能好奇，为啥有那么多工具还要自己写一个。emmm…
比如我求 Ax=b, 增广矩阵, 比如说我想要$[A,I]$ -> $[I,A^{-1}]$, 工具就用不了了
该代码会被我另一篇博客用到 矩阵分析课后题部分题目之计算机辅助(python)

这个代码的作用

• 通过增广矩阵解多元方程组Ax=b， 或者通过化简阶梯型形状判断有无解
• 通过A I, 求得 I $A^{-1}$, 得到矩阵A的逆
• 化简为阶梯型或最简阶梯型
• 计算行列式(如果是方阵的话), 例子见矩阵分析课后题部分题目之计算机辅助(python)chapter 2
• 克拉默求解, 例子见矩阵分析课后题部分题目之计算机辅助(python)chapter 2 q_20
• …

代码

echelon_form.py

"""
支持分数，主要用来计算阶梯型和最简阶梯型, 行列式， 矩阵的逆, 克拉默求解x
Ax = b
Alter time: 2022.9.25 17:19
"""
from fractions import Fraction
import copy


cnt = 0  # 步骤计数器
det_cof = 1  # 经过行变换导致行列式的变化(乘上这个就不变)


# frac.numerator 整数/ frac.denominator分数
# 创造一个支持分数运算的矩阵
def make_matrix(list_matrix):
    global cnt
    global det_cof
    cnt = 0
    det_cof = 1

    fra_matrix = []
    for i in range(len(list_matrix)):
        row = []
        for j in list_matrix[i]:
            if not isinstance(j, Fraction):
                row.append(Fraction(j))
            else:
                row.append(j)
        fra_matrix.append(row)
    show(fra_matrix, 'process {}. former matrix:'.format(cnt))
    return fra_matrix


# 自动增加单位矩阵的 “增广”矩阵, 主要是为了求逆
def make_A_I_matrix(a):
    assert len(a) == len(a[0]), 'len(a) != len(a[0])'
    n = len(a)
    for row in range(n):
        for col in range(n):
            if (row == col):
                a[row].append(1)
            else:
                a[row].append(0)
    return make_matrix(a)


# 展示矩阵
def show(matrix, information_str='no tips', show_procedure=True, cal_det=False):
    global cnt
    if show_procedure:
        print(information_str)
        print('-' * len(matrix[0]) * 10)
        for row in matrix:
            for item in row:
                print('{:>10}'.format(str(item)), end='')
            print()
        if cal_det:
            print('and the cofficient before the det(A) after the row operation is {}|A|'.format(
                det_cof
            ))
        print()
    cnt += 1


# 交换两行 (行数下标从1开始, 即 m * n的矩阵, 行数下标从1 ~ m)
def change_row(matrix, row1_num, row2_num, show_procedure, cal_det=False):
    global cnt
    global det_cof
    matrix[row1_num],  matrix[row2_num] = matrix[row2_num], matrix[row1_num]
    det_cof *= -1    # 交换两行，行列式改变符号
    show(matrix, 'process {}. change row {} and row {}'.format(cnt, row1_num, row2_num), show_procedure, cal_det)


# 化简某一行 (行数下标从0开始, 即 m * n的矩阵, 行数下标[0, m-1] )
def simplify(matrix, row_num, pivot_col, reduce_echelon_form=False, show_procedure=True, cal_det=False):
    global cnt
    global det_cof
    # row_num -= 1
    pivot = matrix[row_num][pivot_col]
    for j in range(len(matrix[0])):
        matrix[row_num][j] /= pivot
    det_cof *= pivot    # 某一行除以多少就要乘回来

    # 化简为 最简行列式
    if reduce_echelon_form:
        for i in range(len(matrix)):
            if i != row_num:
                times = matrix[i][pivot_col]
                for j in range(len(matrix[0])):
                    matrix[i][j] -= matrix[row_num][j] * times
    else: # 只是化为行列式
        for i in range(row_num + 1, len(matrix)):
            times = matrix[i][pivot_col]
            for j in range(len(matrix[0])):
                matrix[i][j] -= matrix[row_num][j] * times
    show(matrix, 'process {}. simplify row {}'.format(cnt, row_num), show_procedure, cal_det)


# 化为阶梯型
def echelon_form(a, cc_num=None, reduce_echelon_form=False, cal_det=False, show_procedure=True):
    m, n = len(a), len(a[0])
    # print('s001: ', m, n, cofficient_cols_num)

    if cc_num is None:
        cc_num = n

    prc_row, prc_col = 0, 0
    while prc_row < m and prc_col < cc_num:
        if a[prc_row][prc_col] == 0:
            zero_flag = True
            for f_row in range(prc_row + 1, m):
                if a[f_row][prc_col] != 0:
                    change_row(a, prc_row, f_row, show_procedure, cal_det)
                    zero_flag = False
                    break
            if zero_flag:
                prc_col += 1
                continue
        if show_procedure:
            print('prc_row, prc_col:', prc_row, prc_col)
        simplify(a, prc_row, prc_col, reduce_echelon_form=reduce_echelon_form,
         show_procedure=show_procedure, cal_det=cal_det)
        prc_row += 1
        prc_col += 1
    
    tips_op = '(reduced)' if reduce_echelon_form else ''

    if prc_row == m and prc_col == cc_num:
        print('The matrix has been changed to {} echelon form'.format(tips_op))
    else:
        print('The matrix is singular')

    if cal_det:
        if len(a) != len(a[0]):
            "you're trying to cal a det of matrix, but this matrix is not square matrix"
            return 'error'
        if prc_row == m and prc_col == cc_num:
            det_result = det_cof
            for i in range(0, m):
                det_result *= a[i][i]
                if det_result == 0:
                    break
            print('the det of this matrix is {}.'.format(det_result))
        else:
            print('Because the matrix is singualr, so its det is 0.')
        return det_result


# 通过Guass-Jordon方法求逆
def cal_inv_by_guass_jordan(a):
    cofficient_cols_num = len(a[0])

    a = make_A_I_matrix(a)
    
    echelon_form(a, cofficient_cols_num, reduce_echelon_form=True,
     cal_det=False, show_procedure=True)


# 将A的某一行进行替换
def replace_A_by_b_at_col_i(a, b, col_i):
    r = copy.deepcopy(a)
    for row in range(len(a)):
        r[row][col_i] = b[row]
    return r


# 使用克拉默计算
def cal_x_by_Cramer(a, b):
    x_lt = []
    cofficient_cols_num = len(a)
    ma = make_matrix(a)
    x_lt.append(echelon_form(ma, cofficient_cols_num, reduce_echelon_form=False,
     cal_det=True, show_procedure=False))
    if x_lt[0] == 'error':
        raise Exception('cal det error')
    
    for col in range(len(a[0])):
        a_im1 = replace_A_by_b_at_col_i(a, b, col)
        ma = make_matrix(a_im1)
        x_lt.append(echelon_form(ma, cofficient_cols_num, reduce_echelon_form=False,
        cal_det=True, show_procedure=False))

    print('\n------ Cramer result ------')
    print('det(A) = {}'.format(x_lt[0]))
    for i in range(1, len(x_lt)):
        print('det(A_{}) = {}, x_{} = ({}) / ({}) = {}'.format(
            i, x_lt[i], i, x_lt[i], x_lt[0], round(float(x_lt[i]) / float(x_lt[0]), 2)))
    print('------ Cramer end ------')


def matrix_mul(a, b, show_result=True):
    r1, c1 = len(a), len(a[0])
    r2, c2 = len(b), len(b[0])
    assert c1 == r2, 'the col of matrix a is not equal to the row of matrix b'

    result = [[0 for col in range(c2)] for row in range(r1)]

    for i in range(r1):
        for j in range(c2):
            for k in range(c1):
                result[i][j] += a[i][k] * b[k][j]
    
    print('-' * c2 * 4)
    for row in range(r1):
        for col in range(c2):
            print('{}'.format(result[row][col]))
    print('-' * c2 * 4)
    return result


if __name__ == '__main__':
    a = [[1, 2, 1, 3], [2, 5, -1, -4], [3, -2, -1, 5]]
    a = make_matrix(a)
    cnt = 0  # 步骤计数器

    cofficient_cols_num = 3

    echelon_form(a, cofficient_cols_num, reduce_echelon_form=True, cal_det=False)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219

例子(chapter 1, question 23 (a)题)输出:

process 0. former matrix:
----------------------------------------
         1         2         1         3
         2         5        -1        -4
         3        -2        -1         5

prc_row, prc_col: 0 0
process 0. simplify row 0
----------------------------------------
         1         2         1         3
         0         1        -3       -10
         0        -8        -4        -4

prc_row, prc_col: 1 1
process 1. simplify row 1
----------------------------------------
         1         0         7        23
         0         1        -3       -10
         0         0       -28       -84

prc_row, prc_col: 2 2
process 2. simplify row 2
----------------------------------------
         1         0         0         2
         0         1         0        -1
         0         0         1         3
The matrix has been changed to (reduced) echelon form
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

如果最后一行代码改一下, 即reduce_echelon_form为False， 我们只想得到阶梯型

echelon_form(a, cofficient_cols_num, reduce_echelon_form=False)
1

输出:

process 0. former matrix:
----------------------------------------
         1         2         1         3
         2         5        -1        -4
         3        -2        -1         5

prc_row, prc_col: 0 0
process 0. simplify row 0
----------------------------------------
         1         2         1         3
         0         1        -3       -10
         0        -8        -4        -4

prc_row, prc_col: 1 1
process 1. simplify row 1
----------------------------------------
         1         2         1         3
         0         1        -3       -10
         0         0       -28       -84

prc_row, prc_col: 2 2
process 2. simplify row 2
----------------------------------------
         1         2         1         3
         0         1        -3       -10
         0         0         1         3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26