2025, Oct 19 18:17

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

Почему порядок переменных в SymPy меняет решения solve для нелинейных систем и как стабилизировать вывод: группировка символов и базисы Грёбнера на примере.

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

Повторяемый пример: та же система, другой порядок переменных

В системе ниже заданы три неизвестных и несколько параметров. Если переставить порядок неизвестных в solve, первые две ветви решения совпадают, а третья отличается по видимой сложности для одной из переменных.

from sympy.solvers import solve as ssolve
import sympy as sp
rho = sp.Symbol('rho')
pf = sp.Symbol('pf')
pb = sp.Symbol('pb')
nb = sp.Symbol('nb')
nt = sp.Symbol('nt')
g = sp.Symbol('g')
pt = sp.Symbol('pt')
delta = sp.Symbol('delta')
kappa = sp.Symbol('kappa')
nt = 1
eq_t = rho*(1 - 2**-nt)*pt + pf*pt*g - pt*(delta + kappa*nt)
eq_b = rho*(1 - 2**-nb)*pb + pf*pb*g - pb*(delta + kappa*nb)
eq_f = rho*pf + pb*rho*2**-nb - g*(pt + pb)*pf - delta*pf
eq_t = sp.simplify(eq_t)
sol = ssolve([eq_b, eq_t, eq_f], [pt, pb, pf], rational=True)
sol2 = ssolve([eq_f, eq_b, eq_t], [pf, pb, pt], rational=True)

Оба вызова возвращают одинаковые две первые ветви: (pf, pb, pt) = (0, 0, 0) и ((-(-2*kappa - 2*delta + rho))/(2*g), 0, (-delta + rho)/g). Третья ветвь также совпадает по pt и pf, но для pb её вид зависит от перестановки: при одном порядке pb раздувается до огромного выражения, при другом остаётся компактным.

В компактной ветви pt = 0, pf = (2**nb*(kappa*nb + delta - rho) + rho)/(2**nb*g), а pb равно (-2**nb*kappa*delta*nb + 2**nb*kappa*nb*rho - 2**nb*delta**2 - 2**nb*rho**2 + 2**(nb + 1)*delta*rho - delta*rho + rho**2)/(2**nb*g*(kappa*nb + delta - rho)).

Что на самом деле происходит

Речь идёт не о различии в математике, а о вопросах упрощения. Оба вида pb эквивалентны, что подтверждается напрямую: упрощение их разности даёт ноль. Такое поведение в SymPy известно; также легко увидеть, что численная подстановка через lambdify при примерных значениях параметров даёт одинаковые результаты, хотя одной числовой проверки недостаточно для доказательства.

Перестановка неизвестных влияет на внутренние шаги исключения и упрощения, из‑за чего итог может выглядеть радикально по‑разному, оставаясь эквивалентным. На самом деле у системы три ветви решения; меняется лишь их запись.

Надёжный способ стабилизировать вывод

Рабочая стратегия — перед решением сгруппировать насыщенные параметрами подвыражения во временные символы. Заменяя повторяющиеся фрагменты короткими обозначениями, вы снижаете разрастание выражений и делаете solve менее чувствительным к перестановкам переменных. После решения просто подставьте обозначения обратно, чтобы вернуть исходные параметры.

from sympy import var, nsimplify
from sympy.solvers import solve as ssolve
from itertools import permutations
eqs = [nsimplify(e, rational=True) for e in [eq_b, eq_t, eq_f]]
K0, K1, K10, K11, K12, K13, K2, K4, K5, K6, K7 = var('K0 K1 K10 K11 K12 K13 K2 K4 K5 K6 K7')
# Сжатые уравнения
eqc = [
    K2*(K6*pb + K7*pb + pb*pf),
    K2*pt*(K4 + pf),
    K13*(K10*pf + K11*pf + K12*pf*(pb + pt) + pb)
]
# Обратная подстановка для плейсхолдеров
M = [
    (K13, K5),
    (K12, -2**nb*g/rho),
    (K11, -2**nb*delta/rho),
    (K10, 2**nb),
    (K7, K1*rho/g),
    (K6, -K0/g),
    (K5, rho/2**nb),
    (K4, -kappa/g - delta/g + rho/g/2),
    (K2, g),
    (K1, 1 - 1/2**nb),
    (K0, kappa*nb + delta)
]
# Решение для всех перестановок становится согласованным
vars3 = (pb, pf, pt)
all_solutions = [ssolve(es, vs) for es in permutations(eqc) for vs in permutations(vars3)]
unique_families = {frozenset(branch) for branches in all_solutions for branch in branches}
# Один показательный прогон и подстановка исходных параметров
sample = all_solutions[0]
sol_expanded = [[comp.simplify() for comp in sp.Tuple(*s).subs(M)] for s in sample]

После обратной подстановки все три ветви получаются компактными и инвариантными к перестановкам — и удовлетворяют исходным уравнениям. Прямая проверка подстановкой с последующим упрощением даёт ноль для всех невязок.

# Проверка: подставляем решения в исходные уравнения
check = {E.subs(dict(zip(vars3, S))).cancel().simplify() for E in eqs for S in sol_expanded}

Набор решений после удаления плейсхолдеров выглядит ровно и стабильно: [0, 0, 0], затем [0, (kappa + delta - rho/2)/g, (-delta + rho)/g] и, наконец, [-(delta - rho)*(2**nb*(kappa*nb + delta) - rho*(2**nb - 1))/(2**nb*g*(kappa*nb + delta - rho)), kappa*nb/g + delta/g - rho/g + rho/(2**nb*g), 0].

Базисы Грёбнера как проверка корректности и сложности

Ещё один взгляд — базисы Грёбнера. Система для groebner посильна, но число операций зависит от перестановки переменных. Для «сжатых» уравнений счётчики такие: [63, 59, 63, 42, 32, 32], а для исходных — [946, 344, 946, 213, 156, 156]. Это согласуется с наблюдением, что группирование повторяющейся структуры параметров снижает сложность и стабилизирует дальнейшую алгебру.

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

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

Выводы

Помните, что solve может выдавать по‑разному выглядящие ответы при разных порядках переменных, даже если решения совпадают. Используйте simplify для разности выражений, чтобы удостовериться в эквивалентности, когда формы не совпадают внешне. Группируйте тяжёлые по параметрам фрагменты в плейсхолдеры до решения — так ветви получаются стабильными и компактными, независимо от перестановок уравнений и переменных. При необходимости подключайте Groebner, чтобы оценить алгебраическую сложность и проверить, что структура системы обрабатывается согласованно.

Статья основана на вопросе с StackOverflow автора Wilco Verweij и ответе от smichr.