当前位置：首页 > article >正文

基于python的线性代数相关计算

article 2024/11/14 18:39:36

使用 Python 的 sympy 模块进行线性代数相关计算, 包括矩阵创建、运算、化简、求秩、行列式计算、逆矩阵求解、伴随矩阵计算以及特征值和特征向量的计算等内容，还包括方程组求解和矩阵正定等性质的判别。

1. 矩阵创建

# 安装sympy模块: pip install sympy
# 导入sympy模块
from sympy import *

# 创建列向量
vec = Matrix([1, 2, 3])  

# 对列向量进行转置
vec_transpose = vec.T  
print(vec_transpose)  
# Matrix([[1, 2, 3]])

# 创建三行三列的一般矩阵
matrix = Matrix([[1, 3, 7], [2, 4, 8], [3, 5, 9]])  
print(matrix)  
# Matrix([[1, 3, 7], [2, 4, 8], [3, 5, 9]])

# 使用Matrix(m,n,list)形式创建两行三列矩阵
matrix_2 = Matrix(2, 3, [1, 2, 3, 4, 5, 6])  
print(matrix_2)  
# Matrix([[1, 2, 3], [4, 5, 6]])

# 创建三阶单位矩阵
eye_matrix = eye(3)  
print(eye_matrix)  
# Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

# 创建三阶零矩阵
zeros_matrix = zeros(3)  
print(zeros_matrix)  
# Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])

# 创建对角矩阵
diag_matrix = diag(1, 2, 3)  
print(diag_matrix)  
# Matrix([[1, 0, 0], [0, 2, 0], [0, 0, 3]])

# 创建四阶全1矩阵
ones_matrix = ones(4)  
print(ones_matrix)  
# Matrix([[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]])

# 创建2行2列符号矩阵
symbols_matrix = Matrix(2, 2, symbols('a(1:3)(1:3)'))  
print(symbols_matrix)  
# Matrix([[a11, a12], [a21, a22]])

# 创建四阶希尔伯特矩阵
hilbert_matrix = Matrix(4, 4, lambda i, j :1/(i + j + 1))  
print(hilbert_matrix)
# Matrix([[1, 1/2, 1/3, 1/4], [1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7]])

2. 矩阵运算

# 生成两个符号矩阵
A = Matrix(3, 3, symbols('a(1:4)(1:4)'))
B = Matrix(3, 3, symbols('b(1:4)(1:4)'))
print(A)
# Matrix([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]])
print(B)
# Matrix([[b11, b12, b13], [b21, b22, b23], [b31, b32, b33]])

# 矩阵加法
add_result = A + B
print(add_result)  
# Matrix([[a11 + b11, a12 + b12, a13 + b13], [a21 + b21, a22 + b22, a23 + b23], [a31 + b31, a32 + b32, a33 + b33]])

# 矩阵减法
sub_result = A - B
print(sub_result)  
# Matrix([[a11 - b11, a12 - b12, a13 - b13], [a21 - b21, a22 - b22, a23 - b23], [a31 - b31, a32 - b32, a33 - b33]])

# 矩阵乘法
mul_result = A * B
print(mul_result)  
# Matrix([[a11*b11 + a12*b21 + a13*b31,
#  a11*b12 + a12*b22 + a13*b32, a11*b13 + a12*b23 + a13*b33],
#  [a21*b11 + a22*b21 + a23*b31, a21*b12 + a22*b22 + a23*b32, 
# a21*b13 + a22*b23 + a23*b33],
#  [a31*b11 + a32*b21 + a33*b31, a31*b12 + a32*b22 + a33*b32,
#  a31*b13 + a32*b23 + a33*b33]])

# 行阶梯形化简
C = Matrix([[8, 6, 3, 4, 5], [4, 5, 6, 7, 8], [7, 8, 9, 10, 11]])
# 返回行最简矩阵
rref_result = C.rref()[0] 
print(rref_result)  
# Matrix([[1, 0, 0, -4, -8], [0, 1, 0, 7, 14], [0, 0, 1, -2, -5]])

# 查看矩阵的秩
rank_result = C.rank()
print(rank_result)  
# 3

# 计算行列式
D = Matrix(4, 4, lambda i, j :1/(i + j + 1))
det_result = D.det()
print(det_result)  
# 1/6048000

# 逆矩阵求解
inv_result = D.inv()
print(inv_result)  
# atrix([[16, -120, 240, -140], [-120, 1200, -2700, 1680], [240, -2700, 6480, -4200], [-140, 1680, -4200, 2800]])

# 伴随矩阵
adjugate_result = D.adjugate()
print(adjugate_result)
# Matrix([[1/378000, -1/50400, 1/25200, -1/43200], [-1/50400, 1/5040, -1/2240, 1/3600], [1/25200, -1/2240, 3/2800, -1/1440], [-1/43200, 1/3600, -1/1440, 1/2160]])

3. 方程组求解

# 解非齐次线性方程组
A1 = Matrix([[1, 1, 1], [1, 2, 3], [1, 1, -1]])
A2 = Matrix([[1, 1, 1], [1, 2, 3], [1, 1, 1]])
A3 = Matrix([[1, 1, 1], [2, 2, 2], [3, 3, 3]])
B = Matrix([1, 2, 3])

# 求解唯一解的方程组
print(linsolve((A1, B)))  
# {(-1, 3, -1)}

# 求解无解的方程组
print(linsolve((A2, B)))  
# EmptySet

# 求解有无穷解的方程组
print(linsolve((A3, B)))
# {(-tau0 - tau1 + 1, tau0, tau1)}

4. 特征值和特征向量计算

# 计算特征值方法一
A = Matrix([[2, 1, 1], [1, -3, 4], [1, 4, -3]])
Lambda = symbols("lambda")
B = det(Lambda * eye(3) - A)
print(factor(B))  # 分解因式查看特征值
# lambda*(lambda - 3)*(lambda + 7)
print(solve(B))  # 求解特征值
# [-7, 0, 3]

# 计算特征值方法二
print(A.eigenvals())  
# {3: 1, -7: 1, 0: 1}

# 计算特征向量
print(A.eigenvects())  
# [(-7, 1, [Matrix([
# [ 0],
# [-1],
# [ 1]])]), (0, 1, [Matrix([
# [-1],
# [ 1],
# [ 1]])]), (3, 1, [Matrix([
# [2],
# [1],
# [1]])])]

5. 矩阵相似对角化判断

在数学和线性代数中，矩阵相似对角化是指将一个方阵转换为一个对角矩阵，同时保持其特征值不变。如果一个矩阵可以通过相似变换变为对角矩阵，那么这个矩阵被称为可对角化的。对于一个方阵 A，如果存在一个可逆矩阵 P 使得：
P^−1 * A * P=D
其中 D 是一个对角矩阵，那么矩阵 A 就是可对角化的，对角矩阵 D 的对角线上的元素是 A 的特征值。

A = Matrix([[2, 1, 1], [1, -3, 4], [1, 4, -3]])
# 判断是否可以相似对角化
print(A.is_diagonalizable())  
# True

# 求解相似对角化的矩阵
P, B = A.diagonalize()
print(P)
# Matrix([[0, -1, 2], [-1, 1, 1], [1, 1, 1]])

print(B)
# Matrix([[-7, 0, 0], [0, 0, 0], [0, 0, 3]])

# # 验证
print(P**-1 * A * P == B)  
# True

6. 正定、负定、半正定、半负定判别

A = Matrix([[1, 2, 2], [2, 1, 1], [2, 2, 2]])

# eigenvalues函数定义如下：
def eigenvalues(matrix):
    Lambda = symbols("lambda")
    B = det(Lambda * eye(3) - matrix)
    print(factor(B)) 
    print(solveset(B))
    
# 方法一：计算特征值判断
print(eigenvalues(A))  
# lambda*(lambda - 5)*(lambda + 1)
# {-1, 0, 5}
# None

# 方法二：使用函数判断
print(A.is_positive_definite)  
# False
print(A.is_negative_definite)  
# False
print(A.is_positive_semidefinite)  
# False
print(A.is_negative_semidefinite)
# False