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

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 (see e.g. the book Ideals, Varieties, and Algorithms by Cox, Little, Oshea), 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) . \]