An application of the resultant

[This article was first published on Saturn Elephant, 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.

I will soon release an update of qspray on CRAN as well as a new package: resultant. This post shows an application of these two packages.

Consider the two algebraic curves f(x,y)=0 and g(x,y)=0 where f(x,y)=y4y3+y22x2y+x4 and g(x,y)=y2x2.

We will derive their intersection points. First of all, let’s plot them.

f <- function(x, y) {
  y^4 - y^3 + y^2 - 2*x^2*y + x^4
}

g <- function(x, y) {
  y - 2*x^2
}

# contour line for f(x,y)=0
x <- seq(-1.2, 1.2, len = 2000)
y <- seq(0, 1, len = 2000)
z <- outer(x, y, f)
crf <- contourLines(x, y, z, levels = 0)
# contour line for g(x,y)=0
x <- seq(-1, 1, len = 2000)
y <- seq(0, 1.5, len = 2000)
z <- outer(x, y, g)
crg <- contourLines(x, y, z, levels = 0)

Theoretically there is only one contour line for both. But for some technical reasons, crf is split into several pieces:

length(crf)
## [1] 16

So we need a helper function to construct the dataframe that we will pass to ggplot2::geom_path:

intercalateNA <- function(xs) {
  if(length(xs) == 1L) {
    xs[[1L]]
  } else {
    c(xs[[1L]], NA, intercalateNA(xs[-1L]))
  }
}
contourData <- function(cr) {
  xs <- lapply(cr, `[[`, "x")
  ys <- lapply(cr, `[[`, "y")
  data.frame("x" = intercalateNA(xs), "y" = intercalateNA(ys))
}

datf <- contourData(crf)
datg <- contourData(crg)

I also plot the intersection points that we will derive later:

datPoints <- data.frame("x" = c(-0.5, 0, 0.5), "y" = c(0.5, 0, 0.5))

library(ggplot2)
ggplot() +
  geom_path(aes(x, y), data = datf, linewidth = 1, color = "blue") +
  geom_path(aes(x, y), data = datg, linewidth = 1, color = "green") +
  geom_point(aes(x, y), data = datPoints, size = 2)

Now we compute the resultant of the two polynomials with respect to x:

# define the two polynomials
library(qspray)
x <- qlone(1)
y <- qlone(2)
P <- f(x, y)
Q <- g(x, y)
# compute their resultant with respect to x
Rx <- resultant::resultant(P, Q, var = 1) # var=1 <=> var="x"
prettyQspray(Rx, vars = "x")
## [1] "16*x^8 - 32*x^7 + 24*x^6 - 8*x^5 + x^4"

We need the roots of the resultant Rx. I use giacR to get them:

library(giacR)
giac <- Giac$new()
command <- sprintf("roots(%s)", prettyQspray(Rx, vars = "x"))
giac$execute(command)
## [1] "[[0,4],[1/2,4]]"
giac$close()
## NULL

Thus there are two roots: 0 and 1/2 (the output of GIAC also provides their multiplicities). Luckily, they are rational, so we can use substituteQspray to replace y with each of these roots in P and Q. We firstly do the substitution y=0:

Px <- substituteQspray(P, c(NA, "0"))
Qx <- substituteQspray(Q, c(NA, "0"))
prettyQspray(Px, "x")
## [1] "x^4"
prettyQspray(Qx, "x")
## [1] "(-2)*x^2"

Clearly, 0 is the unique common root of these two polynomials. One can conclude that (0,0) is an intersection point.

Now we do the substitution y=1/2:

Px <- substituteQspray(P, c(NA, "1/2"))
Qx <- substituteQspray(Q, c(NA, "1/2"))
prettyQspray(Px, "x")
## [1] "x^4 - x^2 + 3/16"
prettyQspray(Qx, "x")
## [1] "1/2 - 2*x^2"

It is easy to see that x=±1/2 are the roots of the second polynomial. And one can check that they are also some roots of the first one. One can conclude that (1/2,1/2) and (1/2,1/2) are some intersection points.

And we’re done.

To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.

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.

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)