How much pizza and how much frozen yogurt? …with Gröbner bases

[This article was first published on R-bloggers on Mikkel Meyer Andersen, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In a recent blog post I tried to get yacas to solve a system of polynomial equations.
Unfortunately it could not do that, so I solved it numerically instead.
Now it is possible – together with many other systems of polynomial equations thanks to
fixing a small error in yacas.
It has now been fixed, also in Ryacas (development version), so
hurry up and update Ryacas to the latest version 0.9.0.9122.

First I will show how to do it directly in Ryacas (and afterwards show what happens behind the curtains).

library(Ryacas)

Solving a system of equations

We define the system from the recent blog post in R:

L <- "x^2 * (y/4) - a*(3*x + 3*y/2 - 45)"
eq1 <- paste0("D(x) ", L) %>% yac_str()
eq2 <- paste0("D(y) ", L) %>% yac_str()
eq3 <- paste0("D(a) ", L) %>% yac_str()
eqs <- c(eq1, eq2, eq3)
eqs
## [1] "(x*y)/2-3*a"      "x^2/4-(3*a)/2"    "45-(3*x+(3*y)/2)"

As known from the theory of Gröbner bases,
instead of finding roots in this entangled system,
we can find the Gröbner bases and find roots for those instead:

eq_str <- paste0("{", paste0(eqs, " == 0", collapse = ", "), "}") %>% 
  y_fn("Solve", "{x, y, a}")
eq_str
## [1] "Solve({(x*y)/2-3*a == 0, x^2/4-(3*a)/2 == 0, 45-(3*x+(3*y)/2) == 0}, {x, y, a})"
sol_all <- eq_str %>% yac_str()
sol_all
## [1] "{{x==0,y==30,a==0},{x==10,y==10,a==50/3}}"

As seen, we get two solutions. The solution \(a = 0\) is not feasible in our case (read the recent blog post) so we just pick the second:

sol_elible_vars <- sol_all %>% y_fn("Nth", 2) %>% yac_str()
sol_elible_vars
## [1] "{x==10,y==10,a==50/3}"

We have now finally arrived at the answer:
\[
(x, y, a) = (10, 10, 50/3) .
\]

It is also possible to postprocess the answer a bit if that is needed:

# Remove variable names if wanted:
sol_elible <- sol_elible_vars %>% y_rmvars() %>% yac_str()
sol_elible
## [1] "{10,10,50/3}"
# And maybe convert to numbers:
sol_elible %>% yac_expr() 
## expression(c(10, 10, 50/3))
# ...maybe even evaluated:
sol_elible %>% yac_expr() %>% eval()
## [1] 10.00000 10.00000 16.66667

Enter Gröbner bases

We could actually have solved the system manually earlier using Gröbner bases.

As known from the theory of Gröbner bases (but I will not start trying to explain it because
I do not know them nearly as well as I should),
instead of finding roots in this entangled system,
we can find the Gröbner bases and find roots for those instead:

fs <- paste0("{", paste0(eqs, collapse = ", "), "}")
fs
## [1] "{(x*y)/2-3*a, x^2/4-(3*a)/2, 45-(3*x+(3*y)/2)}"
gs <- fs %>% y_fn("Groebner") %>% yac_str() %>% as_r()
gs
## [1] "(x*y)/2-3*a"            "x^2/4+((-3)*a)/2"      
## [3] "(-3)*x+((-3)*y)/2+45"   "9*y*a-90*a"            
## [5] "15*y-y^2/2-6*a"         "(-6)*a^2+100*a"        
## [7] "((-15)*y)/2+(-9)*a+225"

As seen the system is now “triangulated”. So we proceed:

yac_str("Solve((-6)*a^2+100*a == 0, a)")
## [1] "{a==0,a==50/3}"

As mentioned, the solution \(a = 0\) is not feasible in our case (read the recent blog post), so we use \(a = 50/3\) and obtain:

yac_str("WithValue(a, 50/3, Solve(9*y*a-90*a == 0, y))")
## [1] "{y==10}"

So \(y = 10\). And we get

yac_str("WithValue({a, y}, {50/3, 10}, Solve((x*y)/2-3*a == 0, x))")
## [1] "{x==10}"

We have now finally arrived at the answer:
\[
(x, y, a) = (10, 10, 50/3) .
\]

To leave a comment for the author, please follow the link and comment on their blog: R-bloggers on Mikkel Meyer Andersen.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)