如何让SymPy强制处理复数化简？矩阵乘法非对角项无法归零问题

问题描述 #

作为SymPy新手，在进行复数矩阵乘法运算时遇到如下问题：
定义的矩阵M和其共轭转置Minv相乘本应得到单位矩阵，但计算M * Minv后，非对角项并未归零（例如得到1/3 -1/3 eps - 1/3 eps_star）。虽然该表达式实际值为0（利用三角函数恒等式可化简为1/3 - 1/3*[2cos(π/3)] = 0），但调用simplify或nsimplify方法无法自动化简为0。

此矩阵用于基变换实现算子对角化，最终需要计算Minv * H * M，当前问题导致无法直接看出结果为对角矩阵。

原始代码 #

矩阵M与Minv定义 #

import sympy as sp

eps = sp.exp(2 * sp.pi * sp.I / 6)

eps_star = sp.conjugate(eps)

M = 1 / sp.sqrt(6) * sp.Matrix([
    [1, 1, 1, 1, 1, 1],
    [1, eps_star, -eps, -1, -eps_star, eps],
    [1, -eps, -eps_star, 1, -eps, -eps_star],
    [1, eps, -eps_star, -1, -eps, eps_star],
    [1, -eps_star, -eps, 1, -eps_star, -eps],
    [1, -1, 1, -1, 1, -1]
]).transpose()

Minv = sp.conjugate(M).transpose()

哈密顿量H矩阵定义 #

alpha, beta = sp.symbols("alpha beta")
H = sp.Matrix([
    [alpha, beta, 0, 0, 0, beta],
    [beta, alpha, beta, 0, 0, 0],
    [0, beta, alpha, beta, 0, 0],
    [0, 0, beta, alpha, beta, 0],
    [0, 0, 0, beta, alpha, beta],
    [beta, 0, 0, 0, beta, alpha],
])

解决方法 #

问题根源在于SymPy默认不会自动利用6次单位根eps的代数/三角函数恒等式（如eps + sp.conjugate(eps) = 2cos(π/3) = 1），因此需要手动指定化简方式，推荐以下两种方案：

方案1：使用三角函数形式定义eps并调用trigsimp #

将eps用三角函数形式定义，再对矩阵乘积的每个元素应用trigsimp（专门处理三角函数化简的方法）：

import sympy as sp

# 改用三角函数形式定义单位根，便于SymPy识别恒等式
eps = sp.cos(2 * sp.pi / 6) + sp.I * sp.sin(2 * sp.pi / 6)
eps_star = sp.conjugate(eps)

# 矩阵定义保持不变
M = 1 / sp.sqrt(6) * sp.Matrix([
    [1, 1, 1, 1, 1, 1],
    [1, eps_star, -eps, -1, -eps_star, eps],
    [1, -eps, -eps_star, 1, -eps, -eps_star],
    [1, eps, -eps_star, -1, -eps, eps_star],
    [1, -eps_star, -eps, 1, -eps_star, -eps],
    [1, -1, 1, -1, 1, -1]
]).transpose()

Minv = sp.conjugate(M).transpose()

# 对矩阵乘积的每个元素应用三角函数化简
MMinv = (M * Minv).applyfunc(lambda x: sp.trigsimp(x))
print(MMinv)  # 输出单位矩阵

# 处理H矩阵的对角化计算
diag_H = (Minv * H * M).applyfunc(lambda x: sp.trigsimp(x))
print(diag_H)  # 输出对角矩阵

方案2：添加单位根约束进行化简 #

如果坚持用指数形式定义eps，可以在化简时显式传入单位根的恒等式约束：

# 指数形式定义eps
eps = sp.exp(2 * sp.pi * sp.I / 6)
eps_star = sp.conjugate(eps)

# 计算非对角项示例
expr = 1/3 - 1/3*eps - 1/3*eps_star
# 添加eps + eps_star = 1的约束进行化简
simplified_expr = sp.simplify(expr, [sp.Eq(eps + eps_star, 1)])
print(simplified_expr)  # 输出0

# 对矩阵整体应用带约束的化简
MMinv = (M * Minv).applyfunc(lambda x: sp.simplify(x, [sp.Eq(eps + eps_star, 1)]))

原理说明 #

• trigsimp会自动识别三角函数的恒等关系，将eps + eps_star转换为2cos(π/3)，进而计算出其值为1，从而将非对角项化简为0。
• 显式添加约束则直接告诉SymPy单位根满足的代数关系，让化简器能够利用这些规则完成归零。

内容的提问来源于stack exchange，提问作者Ger