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 (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) . \]