SymPy求解A@B=0类矩阵方程返回空列表无正确解问题

• 符号变量配置错误：定义a,b等符号时添加了commutative=False参数，该参数适用于矩阵等非交换运算场景，标量运算使用该配置会干扰sympy的求解匹配逻辑。
• 数值精度冲突：通过scipy计算得到的特征值是数值近似解，以此构造的K-w²M矩阵不是严格奇异矩阵，得到的两个方程存在微小的数值误差，并非严格相容，sympy的solve默认做严格相等判断，因此判定为无解返回空列表。
• 类型混用问题：直接使用numpy数值数组和sympy符号矩阵做运算，隐式类型转换会引入额外的精度误差和运算逻辑异常。

方案1：修正sympy求解逻辑 #

修改后的代码如下：

from scipy import linalg
import numpy as np
from sympy import symbols, Matrix, solve, N

M = np.diag([2,2])
K = np.array([[600,-300],
              [-300,300]
              ])

M_inv = linalg.inv(M)
la, v = linalg.eig(M_inv @ K)
# 取特征值实部避免复数警告
wn = np.sqrt(np.real(la))
wn.sort()

def matris(w):
    A = K - w**2*M
    return A

matris_1 = matris(wn[0])
matris_2 = matris(wn[1])

# 去掉不必要的非交换配置
a, b = symbols('a b')
fi_matris = Matrix([a, b])

def mod(matris_n):
    # 先将numpy数组转为sympy矩阵避免类型冲突
    sym_mat = Matrix(matris_n)
    eq1 = sym_mat @ fi_matris
    eq2 = eq1.subs(b, 1)
    # 仅取第一个方程求解，规避数值误差带来的方程不相容问题
    sol = solve(eq2[0], a)
    return [float(N(sol[0], 3)), 1.0]

mod_1 = mod(matris_1)
mod_2 = mod(matris_2)
print(mod_1) # 输出[0.618, 1.0]
print(mod_2) # 输出[-1.618, 1.0]

方案2：直接利用已计算的特征向量（更高效） #

你调用linalg.eig得到的v就是特征向量矩阵，每一列对应一个特征值的特征向量，直接做归一化把第二个分量设为1即可，不需要额外用sympy求解：

from scipy import linalg
import numpy as np

M = np.diag([2,2])
K = np.array([[600,-300],
              [-300,300]
              ])

M_inv = linalg.inv(M)
la, v = linalg.eig(M_inv @ K)
# 特征向量归一化：每列除以第二个元素，固定b=1
v_normalized = v / v[1, :]
mod_1 = np.round(v_normalized[:, 0], 3).tolist()
mod_2 = np.round(v_normalized[:, 1], 3).tolist()
print(mod_1) # 输出[0.618, 1.0]
print(mod_2) # 输出[-1.618, 1.0]

内容的提问来源于stack exchange，提问作者buğra caner