Как избежать зависаний в SymPy при поиске комплексных корней

Показываем, почему solveset в SymPy «виснет» на высоких степенях и как быстро визуализировать и найти комплексные корни: domain coloring, nsolve, Poly.nroots.

Когда вы пытаетесь найти и отобразить комплексные корни выражения большой степени в SymPy, на первый взгляд безобидный вызов может «заморозить» весь процесс. Виновником чаще всего оказывается смешение числовых модулей с символьными и принудительная численная оценка «тяжёлых» символьных корней непосредственно перед построением графика. Ниже — краткий разбор, что именно идёт не так и как это исправить, а также надёжные варианты как для визуализации, так и для численного поиска корней.

Постановка задачи

Следующий фрагмент строит полином высокой степени из параметров, вычисленных через math, запрашивает все комплексные корни символьно через solveset, а затем пытается извлечь действительные и мнимые части для отрисовки. Процесс может зависнуть ещё до появления графика.

from __future__ import division
from sympy import *
import numpy as np
import math
import matplotlib.pyplot as mp
u1, u2, w, phi = symbols('u1 u2 w phi')
step = 0.04 * math.pi
lag = 2.0 * math.pi / step
print(lag)
root_set = solveset(w**math.floor(lag + 0.5) + 1 - w, w, domain=Complexes)
print(root_set)
ang = np.linspace(0, 2 * np.pi, 100)
ux = np.cos(ang)
uy = np.sin(ang)
re_vals = [rv.as_real_imag()[0] for rv in root_set]
im_vals = [rv.as_real_imag()[1] for rv in root_set]
mp.figure(figsize=(6, 6))
mp.scatter(re_vals, im_vals, color='green', marker='o', label='Complex Solutions')
mp.plot(ux, uy, 'b--', alpha=0.5, label='Unit Circle Reference')
mp.xlabel('Real Axis')
mp.ylabel('Imaginary Axis')
mp.title('Solutions on Unit Circle')
mp.gca().set_aspect('equal')
mp.grid(True)
mp.legend()
mp.show()

Что именно идёт не так

Базовое выражение — w**floor(lag + 1/2) + 1 - w, где floor(lag + 1/2) вычисляется в 50. Вызов solveset для него возвращает множество из 50 корней в комплексной области. Многие из них — объекты ComplexRootOf, то есть точные, но невычисленные алгебраические числа. Чтобы построить график, их нужно превратить в числа с плавающей точкой, а это требует численной оценки. В данном случае вычисление даже одного такого корня может занять очень много времени. Короткий эксперимент показывает узкое место:

rt_list = list(root_set)
# вычисление даже одного значения может оказаться непомерно медленным
# rt_list[10].n()

Есть и вторая причина трения. Смешивание math (который производит float) с sympy (который умеет сохранять точность) рано переводит SymPy в численные, а не точные вычисления — и это редко желаемо при построении символьных выражений.

Надёжные исправления

Во‑первых, собирайте выражение только из объектов SymPy. Так сохраняется точная арифметика, а символьные операции остаются согласованными и быстрыми. Во‑вторых, избегайте символьного извлечения корней, если вам нужна лишь визуализация или численные приближения. Ниже — два практичных пути.

Быстрая визуализация без вычисления символьных корней

Окрашивание области (domain coloring) мгновенно показывает нули и полюса на комплексной плоскости без какого‑либо символьного решения. В примере ниже используется SymPy Plotting Backend для отображения модуля и наложения единичной окружности.

from sympy import *
from spb import *
xv, yv, w, ang = symbols('xv yv w ang')
# точные параметры SymPy
step_q = Rational(4, 100) * pi
lag_q = 2 * pi / step_q
pow_deg = floor(lag_q + Rational(1, 2))
poly_f = w**pow_deg + 1 - w
scl = 1.05
graphics(
    domain_coloring(poly_f, (w, -scl - scl*I, scl + scl*I), coloring='k', n=500),
    line_parametric_2d(cos(ang), sin(ang), (ang, 0, 2*pi), rendering_kw={'color': 'b'}, use_cm=False, label='Unit Circle', show_in_legend=False),
    grid=False, aspect='equal', size=(10, 10)
)

В этой амплитудной картине нули видны как чёрные точки. Для данного выражения они распределены вокруг единичной окружности; график показывает их все мгновенно, без обращения к solveset.

Численные корни, когда действительно нужны значения

Если нужны численные приближения всех корней, используйте численный решатель с хорошо выбранными начальными оценками или обратитесь к полиномиальному API. Оба подхода в этом случае работают быстро.

Первый способ — запускать nsolve вдоль единичной окружности и накапливать корни:

from sympy import *
from spb import *
w, th = symbols('w th')
step_q = Rational(4, 100) * pi
lag_q = 2 * pi / step_q
pow_deg = floor(lag_q + Rational(1, 2))
poly_f = w**pow_deg + 1 - w
root_count = 50
offset = (2 * pi / root_count) / 2
found_roots = []
for k in range(root_count):
    theta = 2 * pi / root_count * k + offset
    guess = cos(theta) + I * sin(theta)
    try:
        sol = nsolve([poly_f], [w], [complex(guess)])
        found_roots.append(sol[0, 0])
    except:
        pass
scl = 1.05
th = symbols('th')
graphics(
    domain_coloring(poly_f, (w, -scl - scl*I, scl + scl*I), coloring='k', n=500),
    line_parametric_2d(cos(th), sin(th), (th, 0, 2*pi), rendering_kw={'color': 'b'}, use_cm=False, label='Unit Circle', show_in_legend=False),
    complex_points(*found_roots, rendering_kw={'marker': '.', 'color': 'r'}),
    grid=False, aspect='equal', xlabel='Real', ylabel='Imaginary'
)

Второй — привести выражение к полиному и сразу посчитать численные корни:

from sympy import *
w = symbols('w')
step_q = Rational(4, 100) * pi
lag_q = 2 * pi / step_q
pow_deg = floor(lag_q + Rational(1, 2))
poly_f = w**pow_deg + 1 - w
prec_digits = int(pow_deg)
P = Poly(poly_f, w)
roots_numeric = P.nroots(n=prec_digits)
print(roots_numeric)

Почему это важно

Когда SymPy возвращает объекты ComplexRootOf, преобразование их в числа с плавающей точкой может быть крайне медленным для некоторых многочленов высокой степени. Попытка построить такие корни через извлечение действительных и мнимых частей неявно запускает тяжёлую численную стадию — процесс выглядит как зависание. Создавая выражения из примитивов SymPy и вовремя переходя к численным стратегиям, вы полностью избегаете этой ловушки. Практическое наблюдение из опыта: вычисления могут идти очень быстро, если коэффициент при квадратичном члене близок, но не равен, единице, а корни остаются около единичной окружности.

Выводы

Держите символьный и численный миры отдельно: строите выражения с помощью Rational, pi и floor из SymPy, не смешивая их с math. Не настаивайте на точных формулах для корней, если нужен лишь рисунок или значения с плавающей точкой; для этого подходят domain coloring и численный поиск корней. Если требуется весь набор корней, давайте начальные точки для nsolve вокруг единичной окружности или используйте Poly.nroots — так вы быстро получите численные аппроксимации. Этот подход сохраняет точность там, где это полезно, и включает численные методы там, где они действительно окупаются.