How do I plot a general parametric for conic sections in R (ax^2+bxy+cy^2+dx+ey+f=0)?

48 views Asked by At

Given the general form:

A(x^2)+B(x*y)+C(y^2)+D(x)+E(y)+F=0

How can I go about plotting the upper right (positive x,y) quadrant of these conic sections? I found some code to do it for ellipses using a parametric form, but I'm having trouble figuring out the more general solution.

plot.ellipse <- function (a, b, c, d, e, f, n.points = 1000) {
  ## solve for centre
  A <- matrix(c(a, c / 2, c / 2, b), 2L)
  B <- c(-d / 2, -e / 2)
  mu <- solve(A, B)
  ## generate points on circle
  r <- sqrt(a * mu[1] ^ 2 + b * mu[2] ^ 2 + c * mu[1] * mu[2] - f)
  theta <- seq(0, 2 * pi, length = n.points)
  v <- rbind(r * cos(theta), r * sin(theta))
  ## transform for points on ellipse
  z <- backsolve(chol(A), v) + mu
  ## plot points
  # plot(t(z), type = "l", xlim=c(0,1/d), ylim=c(0,1/e))
  return(t(z))
}

When using this for non-ellipses, I get:

Error in chol.default(A) : 
  the leading minor of order 2 is not positive definite

Haven't really found anything that helps me understand if this just requires a tweak or a completely different approach.

1

There are 1 answers

3
Stéphane Laurent On BEST ANSWER

The nature of the conic section is determined by the sign of the discriminant Delta = B² - 4AC. If Delta > 0, this is an ellipse; if Delta < 0, this is a hyperbola; if Delta = 0, this is a parabola.

For an ellipse or a hyperbola, you can use the PlaneGeometry package.

library(PlaneGeometry)

A <- 2
B <- 3
C <- 4
D <- 1
E <- 0
F <- -1

B*B - 4*A*C # < 0 => it is a hyperbola

H <- HyperbolaFromEquation(
  c(Axx = A, Axy = B, Ayy = C, Bx = D, By = E, C = F)
)
H$plot()

Warning

I don't remember whether the equation is

Axx*x² + Axy*x*y + Ayy*y² + Bx*x + By*y + C

or

Axx*x² + 2*Axy*x*y + Ayy*y² + 2*Bx*x + 2*By*y + C

You'll have to check.