This is not my field of expertise, but I’m curious about whether Julia has some package that can solve my problem…
I have some fraction f(x) of polynomials consisting of a numerator polynomial n(x) and a denominator polynomial d(x) [not a transfer function; the polynomials are “symbolic” in nature…]:
f(x) = \frac{n(x)}{d(x)}
I want to simplify this expression as much as possible, i.e., reduce the order of the numerator and denominator polynomial as much as possible when I have as side constraint that x must be the root of polynomial p(x), i.e., under the constraint that:
p(x) = 0
Questions:
How can I do this in a systematic way? [Is there a theory for this… probably.]
Is there a package in Julia that can do this for me?
Motivating example:
Say I have a constraining polynomial:
Without knowing the answer, I would not be able to figure out this simplification. So I’m curious if there is a method/Julia package that can do this for me.
In the case I consider, yes: all three are real. [The problem describes part of a two-parameter cubic Equation of State in Thermodynamics. x is molar volume per co-volume (the unknown), while u and w are parameters characterizing various EoS… Van der Waals cubic from 1873 has u=w=0. Others may have u=1, w=0, u=1+\sqrt{2}, w=1-\sqrt{2}, etc.]
Not sure if this will satisfy you, but here’s a possible solution (for computing the denominator) using ParamPunPam.jl assuming that we’re looking for an expression for f with 1 as the numerator. EDIT: see below message for a more general solution.
System of equations:
{p\left(x\right)} = 0
{n \cdot r} = {d}
{\frac{1}{f}} = r
Script:
using Nemo, ParamPunPam
R_param, (u, w) = QQ["u", "w"]
R, (r, x) = fraction_field(R_param)["r", "x"] # the order is important
p = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
n = x^2 - 2*x - (u + w)
d = (2*x + u) * (x - 1)^2
F = [p, d - n*r]
paramgb(F)
In the REPL:
julia> using Nemo, ParamPunPam
Welcome to Nemo version 0.43.3
Nemo comes with absolutely no warranty whatsoever
julia> R_param, (u, w) = QQ["u", "w"]
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[u, w])
julia> R, (r, x) = fraction_field(R_param)["r", "x"] # the order is important
(Multivariate polynomial ring in 2 variables over fraction field, AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.FracFieldElem{QQMPolyRingElem}}[r, x])
julia> p = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
x^3 - 3*x^2 + (-3*u - 3*w)*x - u^2 - u*w + w
julia> n = x^2 - 2*x - (u + w)
x^2 - 2*x - u - w
julia> d = (2*x + u) * (x - 1)^2
2*x^3 + (u - 4)*x^2 + (-2*u + 2)*x + u
julia> F = [p, d - n*r]
2-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.FracFieldElem{QQMPolyRingElem}}}:
x^3 - 3*x^2 + (-3*u - 3*w)*x - u^2 - u*w + w
-r*x^2 + 2*r*x + (u + w)*r + 2*x^3 + (u - 4)*x^2 + (-2*u + 2)*x + u
julia> paramgb(F)
✓ # Computing specializations.. Time: 0:00:07
2-element Vector{Any}:
x^3 - 3*x^2 + (-3*u - 3*w)*x - u^2 - u*w + w
r - 3*x - u + 1
Here’s a solution that computes both the numerator and denominator simultaneously. Now there are two new variables, n and d, and f is not necessary any more.
Script:
using Nemo, ParamPunPam
R_param, (x, u) = QQ["x", "u"]
R, (n, d, w) = fraction_field(R_param)["n", "d", "w"] # the order is important
P = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
N = x^2 - 2*x - (u + w)
D = (2*x + u) * (x - 1)^2
F = [P, d*N - D*n]
paramgb(F)
julia> using Nemo, ParamPunPam
Welcome to Nemo version 0.43.3
Nemo comes with absolutely no warranty whatsoever
julia> R_param, (x, u) = QQ["x", "u"]
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, u])
julia> R, (n, d, w) = fraction_field(R_param)["n", "d", "w"] # the order is important
(Multivariate polynomial ring in 3 variables over fraction field, AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.FracFieldElem{QQMPolyRingElem}}[n, d, w])
julia> P = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
(-3*x - u + 1)*w + x^3 - 3*x^2 - 3*x*u - u^2
julia> N = x^2 - 2*x - (u + w)
-w + x^2 - 2*x - u
julia> D = (2*x + u) * (x - 1)^2
2*x^3 + x^2*u - 4*x^2 - 2*x*u + 2*x + u
julia> F = [P, d*N - D*n]
2-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.FracFieldElem{QQMPolyRingElem}}}:
(-3*x - u + 1)*w + x^3 - 3*x^2 - 3*x*u - u^2
(-2*x^3 - x^2*u + 4*x^2 + 2*x*u - 2*x - u)*n - d*w + (x^2 - 2*x - u)*d
julia> paramgb(F)
✓ # Computing specializations.. Time: 0:00:08
2-element Vector{Any}:
w + (-1//3*x^3 + x^2 + x*u + 1//3*u^2)//(x + 1//3*u - 1//3)
n - 1//3//(x + 1//3*u - 1//3)*d
Sidenote: it seems a bit arbitrary that it’s necessary to leave w as a variable instead of a parameter. This is because an equation has to have at least one variable. I wonder if there’s a way around this @sumiya11?
NB:
In case it’s not clear, this should be parsed as \frac{\frac{1}{3}}{x + {\frac{1}{3} \cdot u} - \frac{1}{3}}
FTR, exchanging the u and w above results in a different solution. Thus I guess it always makes sense to try several different possibilities for the parameterization and see which solution is best for the problem.
Please do give other examples, in case they’re interesting.
That should give a different solution. When I suggested values for u and w, I confused the values with another parameterization. With u and w, u=w=0 is van der Waals, u=1, w=0 is RK, and u=2, w=-1 is PR.
Groebner bases can do some simplification, as @nsajko suggested.
The idea is to take the expression to be simplified and the known constraints and compute a Groebner basis of that. A simplified expression will be among the elements of the computed basis, and this simplified expression will be equivalent to the original one modulo constraints.
ParamPunPam.jl computes Groebner bases. A subtle moment is the ordering of variables (reordering variables in @nsajko’s code produces different results). @nsajko, I am not 100% sure at the moment, I think your code is a bit incorrect, see Lemma 1.2 in Redirecting.
I would perhaps first use a more low-level package, Groebner.jl, it can be easier to get right:
Among the discovered relations there is n*x + 1//3*n*u - 1//3*n - 1//3*d. This implies d/n = 3x + u - 1 (dividing the relation by n and multiplying by a constant).
System is constructed in the following way:
P is the constraint
d*N - D*n is the expression d/n = N/D = f to be simplified
t*D - 1 is a technicality :^). See the aforementioned paper. Basically, if our expression has denominator D, then we need this artificial constraint.
This approach is fairly general: it should work with any polynomial or rational p(x) and f(x), and the results obtained this way are canonical, though not unique as a matter of course.
The ordering of variables t, n, d, x, u, w in my code above is “the most important variables go first”. The important thing is that t goes first, and that meta-variables n, d go before the actual variables x, u, w. This strategy should work for other examples as well.
I was thinking of Gröbner basis techniques (since it is a formal algorithm related to polynominals), but I don’t really know the theory…
I took a course involving algebra long time ago, but found it very abstract… The professor defined terms recursively… like a semigroup is… then a group is a semigroup with so and so extra properties. Then a field is a group so and so, and a ring is a field…, and a module is a ring… Or whatever. I lost track after some 3 levels :-o
I have the Ideals, Varieties, and Algorithms book by Cox, Little, and O’Shea somewhere; that is the closest to “readable” I’ve found in this field.
Why do I mention this?
to say that when algorithms mention “ring”, etc., I feel like on thin ice. I can probably follow a recipe, but have a limited idea of how things work.
although fascinated by the Gröbner and similar algoriths, I would need examples to figure out how to use it
I didn’t pursue my idea of considering Gröbner basis ideas because…
I thought that this algorithm is mainly for “multivariable” polynomials, and in my case there is only polynomials in “x”… But perhaps the symbolic parameters “u” and “v” can be considered extra variables??
it has never been clear to me whether the original polynomials entered into the algorith are considered “functions”/polynomials (like my numerator and denominator), or “constraints” (polynomials that are equated to zero). And definitely not clear to me those two interpretations can be mixed
Anyways, I really appreciate your responses, and will study your code when I’m back to my laptop.
Any further clarifying comments are also greatly appreciated on my above mentioned uncertainties…
I will also check the proposals on a related fraction I have.
Let me try to see if I understand what you do in Nemo…
R, (t, n, d, x, u, w) = polynomial_ring(QQ, ["t","n","d","x","u","w"])
Here, polynomial_ring is the constructor. Right?
The first argument of the constructor, QQ, tells Nemo to create polynomials over rational fields, i.e., the coefficients for the polynomial can be rational numbers. Right?
The second argument of the constructor, the vector of strings, tells Nemo what “variable” names you want to operate with. Right?
On the LHS (left hand side), R is a “type constructor``”, e.g., I could in principle do R(2) to change argument “2” from being an integer to being an element in the ring of polynomials?? But in your code, you never use R??
On the LHS, the second argument is a tuple of the names you want to use for elements in the ring (??). Do they have to be the same symbols as used in the RHS vector argument of strings?
In your System vector, d and n are unknown polynomials??? in denominator and numerator, respectively??
The computed vector of 13 elements of Groebner basis elements will, amongst others, contain the simplest solution??? What are the other elements of the 13 element response vector??
You say that to make it work, one needs to add the “equation” t*D-1… I assume this means t*D=1. It is not clear to me why this is needed. I assume that t is the multiplicative inverse of D?
n is the numerator, d is the denominator. Groebner bases only work with polynomials, not with rational functions, so we work with the numerator and denominator explicitly.
Yeah, this seems to be necessary, (although I suppose t*d - 1 would be as correct but more efficient than t*D - 1), even though excluding this equation doesn’t result in a wrong answer in this case. Here’s how I explain the extra equation to myself:
The polynomial d*N - D*n represents the equation {d \cdot N} = {D \cdot n}, which I initially thought would be equivalent to \frac{n}{d} = \frac{N}{D}. However it’s not equivalent, because the former has a symmetry between n and d which the latter lacks: in the latter, the denominator can’t be zero, while the numerator can. The additional equation {t \cdot d} = 1, represented as the polynomial t*d - 1 encodes this additional necessary constraint without constraining d otherwise. Substituting d = 0 into the equation {t \cdot d} = 1 results in the falsity 0 = 1.
On the LHS, the second argument is a tuple of the names you want to use for elements in the ring (??). Do they have to be the same symbols as used in the RHS vector argument of strings?
In LHS, one can put any identifiers – this is variable assignment.
In RHS, the strings in polynomial_ring tell Nemo how to print the variables on the screen. For example,
julia> ring, vars = polynomial_ring(QQ, ["a","b","c"])
julia> vars[1] * vars[2] + vars[3]
a*b + c
In your System vector, d and n are unknown polynomials??? in denominator and numerator, respectively??
You say that to make it work, one needs to add the “equation” t*D-1… I assume this means t*D=1. It is not clear to me why this is needed. I assume that t is the multiplicative inverse of D?
@nsajko accurately pointed out that any solution in n, d to the system
p(x) = 0
n/d = f(x)
is “a simplification” of f(x) modulo p(x).
For solving such (multivariate) systems Groebner bases is a tool. In order to use GB we only need to encode the above system in the language of polynomial rings and ideals:
The equation p(x) = 0 is polynomial, it carries over as is.
The equation n/d = f(x) is rational, we need to rewrite it in terms of polynomials. It is equivalent to n - d f(x) = 0, d \neq 0. Given that f(x) = N(x) / D(x), it is in turn equivalent to two polynomial equations n D(x) - d N(x) = 0, d * t = 1 (following @nsajko’s justification), which are precisely the equations in my System.
I hope that explains a bit why and how the System is formed.
Any element of a Groebner basis is a consequence of the original system p(x) = 0, n/d = f(x). The equations in the resulting basis that depend on t are useless for our purpose of expressing n, d in terms of x, so we discard them.
So I select the equation n*x + 1//3*n*u - 1//3*n - 1//3*d = 0 from the basis because it contains only the variables of interest and because it is simple.
This equation produces a solution n/d to the original system, but there is no telling if it would be the simplest solution for some practical definition of simplicity. It seems to work well in this example. @BLI If you have other examples, can you please share them, I would be interested to see how GB performs there.
Btw, internally, GB algorithms take combinations of input equations with different coefficients and try to reduce the order, very similar to what you have done by hand in the original post.
An advantage of treating “u” and “v” as symbolic variables as I have done in the code is that I get a simplified expression that is valid for almost all numeric values of “u” and “v”.