前言

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

这个代码的作用

• 通过增广矩阵解多元方程组Ax=b， 或者通过化简阶梯型形状判断有无解
• 通过A I, 求得 I A − 1 A^{-1} 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)

例子(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

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

echelon_form(a, cofficient_cols_num, reduce_echelon_form=False)

输出:

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