求解矩阵行列式方程复根遇阻：mpmath失败、sympy卡顿

算力阿峰

2026-7-26

求解矩阵行列式方程的复根问题

我尝试求解由矩阵行列式得到的方程f的复根，编写了如下Python代码：

from sympy import *
from math import sqrt

def Rec_Coef(n, l, Q, rho, i):
    r1 = (1/2)*(1+sqrt(1-4*Q**2))
    r2 = (1/2)*(1-sqrt(1-4*Q**2))
    
    b = (rho*r1**2)/(r1-r2)
    
    if i == 1:
        q = (1/2)*(3+sqrt(9+16*Q**2*(l-1)*(l+2)))
    if i == 2: 
        q = (1/2)*(3-sqrt(9+16*Q**2*(l-1)*(l+2)))
        
    A = l*(l+1)    
        
    alpha = (n**2 + 2*(b+1)*n + 2*b + 1)*r1
    beta = (-2 + r2)*n**2 + (-2 - 2*b*(2-r2) - 4*rho*r1**2 + 6*r2)*n + (q - 2*(rho**2)*(r1**3) - 4*rho*r1**2*(1+b) - 2*b**2*(2-r1**(-1)) - r1*A - (3-2*b)*r2)
    gamma = (1+r2)*n**2 + (2*rho*r1*(1+2*r2)+2*b*(1+r2)-10*r2)*n + (rho*r1*(2 -12*r2 + (rho+2*b)*(1+2*r2)) -1 -q -2*b*(1+3*r2) + b**2*(16 + 8*r2 - (15 - 38*r2 + 26*r2**2)*r1**(-3)) + (A+13)*r2)
    delta = (-n**2 + 2*(3 - rho -b)*n - (9 + 4*rho*b - 6*b - 6*rho))*r2
    
    return alpha, beta, gamma, delta

N=4

α = []
β = []
γ = []
δ = []

l = 2 
Q = 0.2
i = 2

rho = symbols("𝜌")

for j in range(N):
    alpha, beta, gamma, delta = simplify(Rec_Coef(j, 2, 0.2, rho, 2))
    α.append(alpha)
    β.append(beta)
    γ.append(gamma)
    δ.append(delta)  

array = [[0]*(N+1) for row in range(N)]

array[0][0] = β[0]
array[0][1] = α[0]
array[1][0] = γ[1]
array[1][1] = β[1]
array[1][2] = α[1]

for i in range(2,N):
    array[i][i-1] = γ[i] - (array[i-1][i-1]/array[i-1][i-2])*δ[i]
    array[i][i] = β[i] - (array[i-1][i]/array[i-1][i-2])*δ[i]
    array[i][i+1] = α[i] 
    
matrix_aux = simplify(Matrix(array))

M_aux = [[0]*(N) for row in range(N)] 

for i in range(N):
    for j in range(N):
        M_aux[i][j]=(array[i][j])  

matrix = Matrix(M_aux)

f = matrix.det()

我尝试了两种求解复根的方法，但都遇到了问题：

第一种方法：mpmath数值求解

使用mpmath将符号方程转为数值方程，以参考值(-0.17792, -0.74734j)为初始点，代码如下：

import mpmath as mp

f = lambdify(rho, matrix.det())

mp.findroot(f, -0.17-0.75j)

但出现错误：

ValueError: Could not find root within given tolerance. (2.01839854512653723801e-14 > 2.16840434497100886801e-19) Try another starting point or tweak arguments.

第二种方法：sympy符号求解

使用sympy的solve函数求解：

f = matrix.det()
sol = solve(f, rho)
sol = [s.n() for s in sol]
for s in sol:
    print(s)

即使设置N=4，代码仍长时间加载无结果。

解决建议

1. 优化符号计算流程

替换math.sqrt为sympy.sqrt：math.sqrt返回浮点数会引入精度损失，改用sympy的符号平方根，保持表达式精确性。
提前计算常数项：r1、r2、q由固定参数Q=0.2、l=2、i=2计算而来，可在函数外提前计算，避免重复运算：
```
Q = 0.2
l = 2
i = 2
r1 = (1/2)*(1+sqrt(1-4*Q**2))
r2 = (1/2)*(1-sqrt(1-4*Q**2))
if i == 2:
    q = (1/2)*(3-sqrt(9+16*Q**2*(l-1)*(l+2)))
```
再将这些值作为参数传入Rec_Coef函数。
简化矩阵构造：直接构造N×N矩阵，避免多余的matrix_aux和M_aux转换，减少符号运算量。
简化行列式表达式：对f = matrix.det()执行simplify(f)或factor(f)，降低后续求解的复杂度。

2. 改进mpmath数值求解

使用更精确的初始点：直接用参考值-0.17792 - 0.74734j作为初始点，而非近似值。
调整求解参数：增大容差或增加迭代次数，例如：
```
mp.findroot(f, -0.17792 - 0.74734j, maxsteps=1000, tol=1e-12)
```
指定lambdify后端：用sympy的lambdify指定mpmath作为后端，保证数值计算一致性：
```
f = lambdify(rho, simplify(matrix.det()), 'mpmath')
```

3. 使用sympy数值求解替代符号求解

符号求解solve对于高次多项式效率极低，改用数值求解方法：

单根求解：nsolve
直接传入初始点求解指定根：

from sympy import nsolve
sol = nsolve(simplify(f), rho, -0.17792 - 0.74734j)
print(sol.n(10))  # 打印10位精度结果

所有根求解：nroots
将行列式转为多项式后，用nroots数值求解所有根：

poly_f = Poly(simplify(f), rho)
roots = poly_f.nroots()
for root in roots:
    print(root.n(10))

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