给定R中的一般方程，如何绘制椭圆？

Question

椭圆一般方程：

a * x ^ 2 + b * y ^ 2 + c * x * y + d * x + e * y + f = 0

Answer 1

你可以使用我的包PlaneGeometry（希望很快在 CRAN 上）：

library(PlaneGeometry)

ell <- EllipseFromEquation(A = 4, B = 2, C = 3, D = -2, E = 7, F = 1)
box <- ell$boundingbox()
plot(NULL, asp = 1, xlim = box$x, ylim = box$y, xlab = NA, ylab = NA)
draw(ell, col = "yellow", border = "blue", lwd = 2)

Answer 2

我们可以从ellipse的parametric方程开始（以下来自*），我们需要5个参数：中心(xc, yc)或(h,k)另一种表示法，轴长度a, b和x 轴与主轴phi 或tau 之间的角度以另一种表示法表示。

xc <- 1 # center x_c or h
yc <- 2 # y_c or k
a <- 5 # major axis length
b <- 2 # minor axis length
phi <- pi/3 # angle of major axis with x axis phi or tau

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
plot(x,y,pch=19, col='blue')

现在，如果我们想从 cartesian conic 等式开始，这是一个两步过程。

1. 将cartesian方程转换为polar（parametric），形式我们可以使用下面的方程，首先使用下图中的5个方程得到5个参数（取自http://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf ，详细的数学可以在那里找到）。

2. 使用获得的参数绘制椭圆，如上所示。

对于第 (1) 步，我们可以使用以下代码（当我们知道 A,B,C,D,E,F 时）：

M0 <- matrix(c(F,D/2,E/2, D/2, A, B/2, E/2, B/2, C), nrow=3, byrow=TRUE)
M <- matrix(c(A,B/2,B/2,C), nrow=2)
lambda <- eigen(M)$values
abs(lambda - A)
abs(lambda - C) 

# assuming abs(lambda[1] - A) < abs(lambda[1] - C), if not, swap lambda[1] and lambda[2] in the following equations:

a <- sqrt(-det(M0)/(det(M)*lambda[1]))  
b <- sqrt(-det(M0)/(det(M)*lambda[2]))
xc <- (B*E-2*C*D)/(4*A*C-B^2)
yc <- (B*D-2*A*E)/(4*A*C-B^2)
phi <- pi/2 - atan((A-C)/B)*2

对于第 (2) 步，使用以下代码：

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
plot(x,y,pch=19, col='blue')

Answer 3

当您知道椭圆的中心轴和主轴时，另一个答案向您展示了如何绘制椭圆。但它们在一般椭圆方程中并不明显。所以在这里，我将从头开始。

省略数学推导，您需要从以下等式求解中心：

（哎呀：应该是“生成v”而不是“生成u”；我无法修复它，因为原来的LaTeX 现在丢失了，我不想再次输入...... )

这是一个 R 函数：

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")
  }

几点说明：

1. 参数a, b, ..., f 有条件，以确保方程是椭圆而不是其他东西（比如抛物线）。所以，不要传入任意参数值来测试。事实上，从等式中可以大致看出这样的要求。比如矩阵A必须是正定的，所以a > 0和det(A) > 0；还有，r ^ 2 > 0。
2. 我使用了 Cholesky 分解，因为这是我最喜欢的。然而，最漂亮的结果来自特征分解。这部分我不再赘述。如果您对此感兴趣，请阅读我的另一个答案Obtain vertices of the ellipse on an ellipse covariance plot (created by car::ellipse)。有漂亮的图来说明 Cholesky 分解和特征分解的几何。