# Symbolics + Gröbner basis
**URL:** https://discourse.julialang.org/t/symbolics-grobner-basis/127706
**Category:** Modelling & Simulations
**Tags:** polynomials
**Created:** [April 4, 2025, 1:55pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706 "2025-04-04T13:55:27Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 4, 2025, 1:55pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/1 "2025-04-04T13:55:28Z")
I’m trying to find the roots of some multivariable polynomials with parameters. Just to illustrate the problem, suppose I have 2 equations, say:
0 = x^2y + 2xy - a \\ 0 = xy - xy^3
So I try to do the following:
```julia
using Symbolics
using Groebner
@variables x, y
eq1 = x^2*y + 2*x*y - a
eq2 = x*y - x*y^3
sys = [eq1, eq2]
Symbolics.groebner_basis(sys)
```
When I run this code, Symbolics complains that I have not defined `a`.
Fine, I understand that. But if I add `a` in the `@variables` list, how will the algorithm know that I don’t want to solve wrt. a, in other words: that it is `x` and `y` I want to find?
[In my case, I have 20 unknowns and 20 equations… if I add what are really parameters as `@variables`, the function spits out a vector of 168 elements.]
---
### Author:  [@noteas](https://discourse.julialang.org/u/noteas)
#### Post date: [April 4, 2025, 2:12pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/2 "2025-04-04T14:12:25Z")
The theory of Gröbner bases doesn’t have the concept of variables and constants. If you just wanted to find the roots of some multivariable polynomials with parameters, use `symbolic_solve` instead and pass the variables you want to solve (which implicitly utilizes Gröbner bases):
```julia
using Symbolics
using Groebner
@variables x, y, a
eq1 = x^2*y + 2*x*y - a
eq2 = x*y - x*y^3
sys = [eq1, eq2]
symbolic_solve(sys, [x, y])
```
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [April 4, 2025, 2:26pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/3 "2025-04-04T14:26:01Z")
> [@BLI](#):
>
> how will the algorithm know that I don’t want to solve wrt. a, in other words: that it is `x` and `y` I want to find?
I don’t know the specific answer, but as as far as I understand, the answer is basically _by specifying a **monomial ordering** explicitly_. See Groebner.jl. You may need to use it directly, instead of Symbolics.jl, not sure.
There’s also another package that depends on Groebner.jl, ParamPunPam.jl. It aims to support specifying parameters, however as far as I understand it assumes the value of a parameter is never algebraic.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 4, 2025, 2:41pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/4 "2025-04-04T14:41:42Z")
Hm… with my 20 equations and 20 variables, `symbolic_solve` is now into its 18 minutes to find the solution. I use a laptop with the top 2024 AMD processor with 12 cores (don’t know if `symbolic_solve` can utilize the cores, though).
- How large systems is it possible to solve? Is there a limit?
- What is a realistic computational time?
---
### Author:  [@noteas](https://discourse.julialang.org/u/noteas)
#### Post date: [April 4, 2025, 2:48pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/5 "2025-04-04T14:48:01Z")
Even your example took about 1 minute to solve on my 10-year-old computer. I’m not sure how long it’s going to be, but I know that the time complexity of computing Gröbner bases is very very high (Buchberger’s algorithm), although there are faster algorithms like F5 algorithm.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [April 4, 2025, 3:34pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/6 "2025-04-04T15:34:22Z")
The run time also greatly depends on the chosen monomial ordering.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 5, 2025, 7:18am UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/7 "2025-04-05T07:18:27Z")
OK – I had to stop the `symbolic_solve()` algorithm… it had spent more than 12.5 _hours_ without finding the solution during the night.
One pass of running `Symbolics.groebner_basis()` takes less than 7 seconds on my laptop for the 20 equations with 20 unknowns + 10 parameters.
My “toy” example of 2 unknowns and 1 parameter took 13 seconds to solve on my laptop using `symbolic_solve`.
* * *
In my case, I am ok with the solution consisting in a sequence of monovariable polynomials that I can solve using a root-finding algorith, and then substituting in the solutions into the next polynomial.
Perhaps `symbolic_solve()` attempts to actually also solve the polynomials and find explicit solutions?
* * *
As mentioned, the `Symbolics.groebner_basis()` algorithm spits back 168 bases within 7 seconds. I will try to figure out what monomial ordering is.
---
### Author:  [@sumiya11](https://discourse.julialang.org/u/sumiya11)
#### Post date: [April 5, 2025, 11:54am UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/8 "2025-04-05T11:54:35Z")
Hi @BLI ,
To complement a bit on previous comments:
> How large systems is it possible to solve? Is there a limit?
This depends on the structure of the system. One rule of thumb is the Bézout bound.
With 20 equations in 20 variables, say each equation is of total degree 2, there are at most 2^20 solutions (d^n). The solving algorithm needs to compute and represent those in some way.
For example, solving random systems in 30 variables of degree 2 is hopeless with Groebner bases.
> What is a realistic computational time?
Roughly speaking, incrementing the number of variables by one doubles the running time and memory.
> One pass of running `Symbolics.groebner_basis()` takes less than 7 seconds on my laptop for the 20 equations with 20 unknowns + 10 parameters.
This could mean that this system has some nice structure, and the rule I gave above turns out too pessimistic.
> if I add what are really parameters as `@variables` , the function spits out a vector of 168 elements.
> Perhaps `symbolic_solve()` attempts to actually also solve the polynomials and find explicit solutions?
Yes, I think so. Usually, a GB is not enough for getting a useful representation of solutions (i.e., there may be extra steps to pass from a GB with 168 elements to solutions in some form).
> OK – I had to stop the `symbolic_solve()` algorithm… it had spent more than 12.5 _hours_ without finding the solution during the night.
Since the general solver seems stuck, I would perhaps proceed by analysing the system manually: checking the number of solutions and the dimension of the variety from the GB, obtaining useful relations from the GB to simplify the input system, etc.
Feel free to send the system to me ([asdemin\_2@edu.hse.ru](mailto:asdemin_2@edu.hse.ru)) or post it here, I can have a look.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 5, 2025, 12:41pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/9 "2025-04-05T12:41:10Z")
OK – I’ll post it here. It is not a revolutionary problem…
```julia
using Symbolics
using Groebner
# Variables, i.e., unknowns
@variables V_L, V_g, n_H2O_L, n_H2O_g, n_1_L, n_2_L, n_1_g, n_2_g, n_L, n_g
@variables p, p_H2O, p_1, p_2, y_H2O, y_1, y_2, x_H2O, x_1, x_2
#
# Parameters, i.e., known values
@variables n_H2O, n_1, n_2, V, p_H2O_sat, H_1, H_2, Vt_L, R, T
eq1 = n_H2O_L + n_H2O_g - n_H2O
eq2 = n_1_L + n_1_g - n_1
eq3 = n_2_L + n_2_g - n_2
eq4 = n_H2O_g + n_1_g + n_2_g - n_g
eq5 = n_H2O_L + n_1_L + n_2_L - n_L
eq6 = p_H2O + p_1 + p_2 - p
eq7 = p_H2O*V_g - n_H2O_g*R*T
eq8 = p_1*V_g - n_1_g*R*T
eq9 = p_2*V_g - n_2_g*R*T
eq10 = x_H2O*n_L - n_H2O_L
eq11 = x_1*n_L - n_1_L
eq12 = x_2*n_L - n_2_L
eq13 = y_H2O*n_g - n_H2O_g
eq14 = y_1*n_g - n_1_g
eq15 = y_2*n_g - n_2_g
eq16 = p_H2O - x_H2O*p_H2O_sat
eq17 = p_1 - x_1*H_1
eq18 = p_2 - x_2*H_2
eq19 = V_L - n_H2O_L*Vt_L
eq20 = V_L + V_g - V
eqs_VT_flash = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20]
vars_VT_flash = [V_L, V_g, n_H2O_L, n_H2O_g, n_1_L, n_2_L, n_1_g, n_2_g, n_L, n_g, p, p_H2O, p_1, p_2, y_H2O, y_1, y_2, x_H2O, x_1, x_2]
# sol_VT_flash = symbolic_solve(eqs_VT_flash, vars_VT_flash)
sol_VT_flash = Symbolics.groebner_basis(eqs_VT_flash)
```
Essentially, it is an equilibrium thermodynamics problem, a so-called VT Flash problem.
The “parameters” are known values. What I am looking for is the solution of the unknowns. If I could find what I call n\_H2O\_L (amount of water in the liquid phase) and n\_1\_g, n\_2\_g (amount of components 1 and 2 in the gas phase), it is trivial to find the rest.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 5, 2025, 12:49pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/10 "2025-04-05T12:49:54Z")
This should end up down the symbolic linear solver path and not hit the Groebner stuff at all? That’s the real bug here.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 5, 2025, 2:09pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/11 "2025-04-05T14:09:23Z")
I’m not claiming a bug. Perhaps it ends up with a 5th order polynomial or higher that cannot be solved to an explicit solution. But it would be great if `solve_symbolic` could return something in less time. I wouldn’t mind if it was in the form of a monovariable polynomial that needs to be solved numerically (i.e., an implicit expression) + additional polynomials that become monovariable by successively inserting the solution to a prior polynomial.
If it is a bug, or a short-coming, great if my example can inspire further development.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 5, 2025, 2:15pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/12 "2025-04-05T14:15:08Z")
I’m claiming a bug. It should prove it’s linear and short circuit to the symbolic linear solver which is much faster.
---
### Author:  [@sumiya11](https://discourse.julialang.org/u/sumiya11)
#### Post date: [April 5, 2025, 2:36pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/13 "2025-04-05T14:36:26Z")
> [@ChrisRackauckas](#):
>
> It should prove it’s linear and short circuit to the symbolic linear solver which is much faster.
Is it linear? For example, `eq7` is quadratic in `p_H2O`, `V_g`, which are `@variables`.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 5, 2025, 2:43pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/14 "2025-04-05T14:43:34Z")
> [@sumiya11](#):
>
> p\_H2O
Oh I thought that was one of the parameters not differentiated 😅. The naming scheme is a bit rough.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 5, 2025, 3:35pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/15 "2025-04-05T15:35:21Z")
With some simplifying assumptions, it can be made linear. But I am trying to avoid those assumptions.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 5, 2025, 3:43pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/16 "2025-04-05T15:43:47Z")
> [@ChrisRackauckas](#):
>
> The naming scheme is a bit rough.
Sorry about that 🙂 . `n` is amount (mol), `p` is pressure, `V` is volume, `x` is liquid phase mole fraction, `y` is gas phase mole fraction. Subscripts are `H2O` (water), `1`, `2` (two different light gases). “Superscripts” (the second `_` are `L` for liquid, and `g` for gas. `T` is absolute temperature, `R` is gas constant, etc.
What is known at the outset in `VT` (volume-temperature) “flash” is total volume `V`, temperature `T`, and total amounts of the various components (`n_H2O`, `n_1`, `n_2`). I want to know how the components split in the liquid and gas phase.
Finally, `H_1`, `H_2` (Henry’s constants), and `p_H2O_sat` (water saturation pressure), and `Vt_t_L` (liquid molar volume) are known functions of temperature.
---
### Author:  [@sumiya11](https://discourse.julialang.org/u/sumiya11)
#### Post date: [April 5, 2025, 7:32pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/17 "2025-04-05T19:32:18Z")
@BLI thanks for sharing the equations.
Here is how I would get a univariate polynomial in `n_H2O_L` that depends on parameters:
```julia
using Nemo, Groebner
# 1. Define parameters and variables
_r, (n_H2O, n_1, n_2, V, p_H2O_sat, H_1, H_2, Vt_L, R, T) = Nemo.QQ[:n_H2O, :n_1, :n_2, :V, :p_H2O_sat, :H_1, :H_2, :Vt_L, :R, :T]
_R, (V_L, V_g, n_H2O_L, n_H2O_g, n_1_L, n_2_L, n_1_g, n_2_g, n_L, n_g, p, p_H2O, p_1, p_2, y_H2O, y_1, y_2, x_H2O, x_1, x_2) = fraction_field(_r)[:V_L, :V_g, :n_H2O_L, :n_H2O_g, :n_1_L, :n_2_L, :n_1_g, :n_2_g, :n_L, :n_g, :p, :p_H2O, :p_1, :p_2, :y_H2O, :y_1, :y_2, :x_H2O, :x_1, :x_2]
# 2. Define equations
eq1 = n_H2O_L + n_H2O_g - n_H2O
eq2 = n_1_L + n_1_g - n_1
eq3 = n_2_L + n_2_g - n_2
eq4 = n_H2O_g + n_1_g + n_2_g - n_g
eq5 = n_H2O_L + n_1_L + n_2_L - n_L
eq6 = p_H2O + p_1 + p_2 - p
eq7 = p_H2O*V_g - n_H2O_g*R*T
eq8 = p_1*V_g - n_1_g*R*T
eq9 = p_2*V_g - n_2_g*R*T
eq10 = x_H2O*n_L - n_H2O_L
eq11 = x_1*n_L - n_1_L
eq12 = x_2*n_L - n_2_L
eq13 = y_H2O*n_g - n_H2O_g
eq14 = y_1*n_g - n_1_g
eq15 = y_2*n_g - n_2_g
eq16 = p_H2O - x_H2O*p_H2O_sat
eq17 = p_1 - x_1*H_1
eq18 = p_2 - x_2*H_2
eq19 = V_L - n_H2O_L*Vt_L
eq20 = V_L + V_g - V
eqs_VT_flash = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20]
vars_VT_flash = [V_L, V_g, n_H2O_L, n_H2O_g, n_1_L, n_2_L, n_1_g, n_2_g, n_L, n_g, p, p_H2O, p_1, p_2, y_H2O, y_1, y_2, x_H2O, x_1, x_2]
# 3. Make these variables the smallest
interesting_vars = [n_H2O_L]
# 4. Compute a GB, allowing division by parameters
using ParamPunPam
gb = paramgb(
eqs_VT_flash,
# vars > interesting_vars > params
ordering=DegRevLex(setdiff(vars_VT_flash, interesting_vars)) * DegRevLex(interesting_vars)
)
# 5. Inspect. Note: gb[1] is univariate
map(vars, gb)
#= Output
20-element Vector{Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.FracFieldElem{QQMPolyRingElem}}}}:
[n_H2O_L]
[n_H2O_L, x_2]
[n_H2O_L, x_1]
[n_H2O_L, x_H2O]
[n_H2O_L, y_2]
[n_H2O_L, y_1]
[n_H2O_L, y_H2O]
[n_H2O_L, p_2]
[n_H2O_L, p_1]
[n_H2O_L, p_H2O]
[n_H2O_L, p]
[n_H2O_L, n_g]
[n_H2O_L, n_L]
[n_H2O_L, n_2_g]
[n_H2O_L, n_1_g]
[n_H2O_L, n_2_L]
[n_H2O_L, n_1_L]
[n_H2O_L, n_H2O_g]
[V_g, n_H2O_L]
[V_L, n_H2O_L]
=#
gb[1]
#= Output
n_H2O_L^4 + (n_H2O*p_H2O_sat^2*H_1*Vt_L + n_H2O*p_H2O_sat^2*H_2*Vt_L + n_H2O*p_H2O_sat^2*R*T - 2*n_H2O*p_H2O_sat*H_1*H_2*Vt_L - 2*n_H2O*p_H2O_sat*H_1*R*T - 2*n_H2O*p_H2O_sat*H_2*R*T + 3*n_H2O*H_1*H_2*R*T - n_1*p_H2O_sat^2*R*T + n_1*p_H2O_sat*H_2*R*T - n_2*p_H2O_sat^2*R*T + n_2*p_H2O_sat*H_1*R*T - V*p_H2O_sat^3 + V*p_H2O_sat^2*H_1 + V*p_H2O_sat^2*H_2 - V*p_H2O_sat*H_1*H_2)//(p_H2O_sat^3*Vt_L - p_H2O_sat^2*H_1*Vt_L - p_H2O_sat^2*H_2*Vt_L - p_H2O_sat^2*R*T + p_H2O_sat*H_1*H_2*Vt_L + p_H2O_sat*H_1*R*T + p_H2O_sat*H_2*R*T - H_1*H_2*R*T)*n_H2O_L^3 + (n_H2O^2*p_H2O_sat*H_1*H_2*Vt_L + n_H2O^2*p_H2O_sat*H_1*R*T + n_H2O^2*p_H2O_sat*H_2*R*T - 3*n_H2O^2*H_1*H_2*R*T + n_H2O*n_1*p_H2O_sat^2*R*T - 2*n_H2O*n_1*p_H2O_sat*H_2*R*T + n_H2O*n_2*p_H2O_sat^2*R*T - 2*n_H2O*n_2*p_H2O_sat*H_1*R*T - n_H2O*V*p_H2O_sat^2*H_1 - n_H2O*V*p_H2O_sat^2*H_2 + 2*n_H2O*V*p_H2O_sat*H_1*H_2)//(p_H2O_sat^3*Vt_L - p_H2O_sat^2*H_1*Vt_L - p_H2O_sat^2*H_2*Vt_L - p_H2O_sat^2*R*T + p_H2O_sat*H_1*H_2*Vt_L + p_H2O_sat*H_1*R*T + p_H2O_sat*H_2*R*T - H_1*H_2*R*T)*n_H2O_L^2 + (n_H2O^3*H_1*H_2*R*T + n_H2O^2*n_1*p_H2O_sat*H_2*R*T + n_H2O^2*n_2*p_H2O_sat*H_1*R*T - n_H2O^2*V*p_H2O_sat*H_1*H_2)//(p_H2O_sat^3*Vt_L - p_H2O_sat^2*H_1*Vt_L - p_H2O_sat^2*H_2*Vt_L - p_H2O_sat^2*R*T + p_H2O_sat*H_1*H_2*Vt_L + p_H2O_sat*H_1*R*T + p_H2O_sat*H_2*R*T - H_1*H_2*R*T)*n_H2O_L
=#
```
Some comments:
- The ordering of indeterminates is: vars \> interesting vars \> params, + allowing division by params. This can be modeled differently, but I guess this approach works here.
- In 3., mark several variables as interesting to uncover relations between them (for example, [n\_H2O\_L, n\_1\_g, n\_2\_g]).
- I use `ParamPunPam.paramgb` for computing GB, which @nsajko mentioned.
- In 4., swap `paramgb` for `groebner` if guaranteed result is needed (`paramgb` is probabilistic).
- There are 4 solutions, which agrees with the degree of univariate polynomial:
```julia
julia> length(quotient_basis(map(identity, gb)))
4
```
- From the univariate polynomial, one solution is `n_H2O_L = 0` (I wonder if this has physical meaning).
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [April 5, 2025, 7:49pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/18 "2025-04-05T19:49:14Z")
> [@sumiya11](#):
>
> ```julia
> # 5. Inspect. Note: gb[1] is univariate
> map(vars, gb)
>
> ```
I don’t think the value of `vars` is defined?
---
### Author:  [@sumiya11](https://discourse.julialang.org/u/sumiya11)
#### Post date: [April 5, 2025, 7:50pm UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/19 "2025-04-05T19:50:15Z")
It comes from Nemo.jl.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [April 6, 2025, 7:24am UTC](https://discourse.julialang.org/t/symbolics-grobner-basis/127706/20 "2025-04-06T07:24:06Z")
Thanks a lot!!
I’ll need some time to digest it. Seems like a 4th order polynomial in `n_H2O_L`, right? I’ll have to check how the other variables relate to `n_H2O_L`.
Of course, a fourth order polynomial admits an explicit solution, so I guess in principle, `symbolic_solve` could find the solution. But perhaps it is necessary to specify lex ordering, etc. to find the solution.
I tried to scale down the problem to water + one other gas and toss it to `symbolic_solve`, but still no solution after 7.5 hours. In real life problems involving water electrolysis (as an example), there could probably be water + 4 gases (I included water + 2 gases): hydrogen, oxygen, nitrogen, argon. Possibly more. At some stage, it would probably be better to simply solve it as a numeric problem with iteration on implicit equations.
[Next page](https://discourse.julialang.org/t/symbolics-grobner-basis/127706.md?page=2)