[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) = y^4 – y^3 + y^2 – 2x^2y + x^4 \quad \text{ and } \quad g(x, y) = y – 2x^2.$ 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 $$R_x$$. 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= \pm 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.