SymPy求解特定高次复多项式根时耗时过长的异常情况求助

我刚接触SymPy和Python，本来应该用R的，但想试试AI生成的代码。其他测试多项式运行都没问题，但这次在生成符号结果后，绘图前好像卡住了。下面是我的代码，注释掉的那些多项式几分钟就能跑完，有没有大佬能给点思路？谢谢！

初始问题代码 #

from __future__ import division
from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
init_printing()
import numpy as np
import math
import matplotlib.pyplot as plt

delta=.04*math.pi
delay=2.0*math.pi/delta
print ( delay )

# 以下注释的多项式可正常运行
#complex_solutions = solveset(z**2-1-2*z+delta*delta*(z**delay), z, domain=Complexes)
#complex_solutions = solveset(z**2-1-2*z+delta*delta*(z**math.floor(delay+.5)), z, domain=Complexes)
complex_solutions = solveset(z**math.floor(delay+.5)+1-z, z, domain=Complexes)
#complex_solutions = solveset(z**2-1-2*z+z**10, z, domain=Complexes)
print ( complex_solutions )

t = np.linspace(0, 2 * np.pi, 100) # Angle from 0 to 2pi
x_circle = np.cos(t)
y_circle = np.sin(t)

real_parts = [sol.as_real_imag()[0] for sol in complex_solutions]
imag_parts = [sol.as_real_imag()[1] for sol in complex_solutions]

plt.figure(figsize=(6,6))
plt.scatter(real_parts, imag_parts, color='green', marker='o', label='Complex Solutions')
plt.plot(x_circle, y_circle, 'b--', alpha=0.5, label='Unit Circle Reference')
plt.xlabel('Real Axis')
plt.ylabel('Imaginary Axis')
plt.title('Solutions on Unit Circle')
plt.gca().set_aspect('equal')
plt.grid(True)
plt.legend()
plt.show()

问题分析 #

后来我琢磨明白卡住的原因了：当多项式次数极高时（比如这里math.floor(delay+.5)算出来是157次），solveset尝试符号求解高次多项式的复根会异常耗时。根据阿贝尔定理，次数≥5的多项式没有通用根式解，SymPy的符号求解引擎会在这种情况下做大量无效尝试，直接导致程序卡死。

亲测有效的解决方法 #

我换了两种数值求解的思路，都能快速出结果，完美绕开符号求解的坑：

方法1：用Poly.nroots()做多项式数值根求解

这是SymPy专门为多项式优化的数值根求解函数，处理高次多项式的效率拉满：

# 承接前面的代码，expr为你的高次多项式
from sympy import Poly
nlead = math.floor(delay+.5)
expr = z**nlead + 1 - z
p = Poly(expr, z)
complex_solutions = p.nroots(n=nlead)  # n参数指定要求解的根的数量
print(complex_solutions)

方法2：用nsolve循环取单位圆初始猜测值

这类多项式的根大多分布在单位圆附近，我们可以在单位圆上均匀取点作为初始猜测，用nsolve逐个求解：

# 承接前面的代码，expr为你的高次多项式
n = nlead  # 多项式次数
solutions = []
for i in range(n):
    # 在单位圆上取初始猜测点，稍微偏移避免落在坐标轴上
    t_angle = 2 * np.pi / n * i + (2 * np.pi / n)/2
    init_guess = cos(t_angle) + I * sin(t_angle)
    try:
        sol = nsolve([expr], [z], [complex(init_guess)])
        solutions.append(sol[0, 0])
        print(i, "\t", sol[0, 0])
    except:
        pass
complex_solutions = solutions

完整修改后可运行代码 #

from __future__ import division
from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
init_printing()
import numpy as np
import math
import matplotlib.pyplot as plt

delta=Rational(4,100)*pi
delay=2/delta*pi
print ( delay )
print ( delay-math.floor(delay) )

nlead=math.floor(delay+.5)
expr=z**nlead +1 -z

# 选择其中一种求解方法即可
# 方法1：Poly.nroots快速求解
from sympy import Poly
p = Poly(expr, z)
complex_solutions = p.nroots(n=nlead)
print(complex_solutions)

# 绘图部分和原代码一致
t = np.linspace(0, 2 * np.pi, 100)
x_circle = np.cos(t)
y_circle = np.sin(t)

real_parts = [sol.as_real_imag()[0] for sol in complex_solutions]
imag_parts = [sol.as_real_imag()[1] for sol in complex_solutions]

plt.figure(figsize=(6,6))
plt.scatter(real_parts, imag_parts, color='green', marker='o', label='Complex Solutions')
plt.plot(x_circle, y_circle, 'b--', alpha=0.5, label='Unit Circle Reference')
plt.xlabel('Real Axis')
plt.ylabel('Imaginary Axis')
plt.title('Solutions on Unit Circle')
plt.gca().set_aspect('equal')
plt.grid(True)
plt.legend()
plt.show()

核心思路就是：高次多项式别死磕符号求解，直接上数值方法就对了！

内容来源于stack exchange