Stabilizing SymPy nonlinear system solutions: variable order effects, placeholders, and Groebner

Learn why SymPy solve yields different forms under variable order, and how placeholders and Groebner bases tame expression swell for stable nonlinear solutions.

When solving nonlinear systems in SymPy, simply reordering the unknowns in solve can sometimes change the form of the result. The solutions remain mathematically equivalent, but one permutation can explode into a huge expression while another stays compact. This guide shows the behavior on a concrete system, explains why the discrepancy shows up in practice, and demonstrates a robust way to stabilize the output using temporary symbol grouping and Groebner-driven structure.

Repro case: same system, different variable order

The system below uses three unknowns and several parameters. Swapping the order of the unknowns in solve yields matching first and second solutions, but the third one differs in apparent complexity for one variable.

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)

Both calls return the same two first solution branches, namely (pf, pb, pt) = (0, 0, 0) and ((-(-2*kappa - 2*delta + rho))/(2*g), 0, (-delta + rho)/g). The third branch also matches for pt and pf but looks different for pb depending on the permutation. In one ordering pb expands into a very large expression, while in the other it appears in a compact form.

For the compact branch, pt = 0 and pf = (2**nb*(kappa*nb + delta - rho) + rho)/(2**nb*g), and pb equals (-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)).

What’s actually happening

The discrepancy is a simplification issue rather than a mathematical difference. The two forms of pb are equivalent, and a direct equivalence check confirms it: simplifying their difference gives zero. This is a known behavior in SymPy, and it can also be seen that numerical evaluation via lambdify with sample parameter values makes the two expressions agree, though numeric checks alone are not proofs.

Changing the order of unknowns affects internal elimination and simplification paths, which in turn can lead to dramatically different but equivalent output expressions. The system really has three solution branches here; the forms just vary with ordering.

A stable way to tame the output

An effective strategy is to group parameter-heavy subexpressions into temporary symbols before solving. By replacing recurring fragments with compact placeholders, you reduce expression swell and make solve less sensitive to variable permutations. After solving, substitute the placeholders back to recover the original parameters.

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')

# Condensed equations
eqc = [
    K2*(K6*pb + K7*pb + pb*pf),
    K2*pt*(K4 + pf),
    K13*(K10*pf + K11*pf + K12*pf*(pb + pt) + pb)
]

# Back-substitution map for placeholders
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)
]

# Solving across all permutations becomes consistent
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}

# One representative solve and substitution back to original parameters
sample = all_solutions[0]
sol_expanded = [[comp.simplify() for comp in sp.Tuple(*s).subs(M)] for s in sample]

After substitution, the three branches are compact and permutation-invariant, and they satisfy the original equations. A direct verification by substitution and simplification yields zero for all residuals.

# Verify solutions solve the original equations
check = {E.subs(dict(zip(vars3, S))).cancel().simplify() for E in eqs for S in sol_expanded}

The solution set, with placeholders removed, appears in a clean, stable form: [0, 0, 0], then [0, (kappa + delta - rho/2)/g, (-delta + rho)/g], and finally [-(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 as a sanity and complexity check

Another angle is to look at Groebner bases. The system is tractable with groebner, but the operation count changes with variable permutations. On the condensed equations the counts are [63, 59, 63, 42, 32, 32], whereas on the raw equations they are [946, 344, 946, 213, 156, 156]. This is consistent with the observation that grouping recurring parameter structure reduces complexity and stabilizes downstream algebra.

Why this matters

In downstream workflows such as computing Jacobians and eigenvalues, having compact, permutation-stable expressions is more than cosmetic. It makes symbolic differentiation faster and less fragile, improves readability for colleagues, and cuts down the chance of misinterpretation caused by expression swell. The condensed approach leads to smaller expressions without altering the mathematics, and equivalence checks via simplify confirm correctness.

Takeaways

Be aware that solve can yield different-looking outputs across variable orderings even when the solutions are the same. Use simplify on differences to certify equivalence when forms look mismatched. Grouping parameter-heavy fragments into placeholders before solving leads to stable, compact branches that are invariant to equation and variable permutations. If needed, use Groebner to gauge algebraic complexity and confirm that the system structure is being handled consistently.