R语言contourLines绘制对称隐函数曲线仅显示一半的问题及解决

要绘制由隐函数 $f(x, y) = y^4 - y^3 + y^2 - 2x^2y + x^4 = 0$ 定义的曲线，该曲线关于$x$轴对称（因为$f(-x, y) = f(x, y)$），但使用contourLines时仅能获取曲线的一半，代码如下：

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

x <- seq(-3, 3, len = 1000)
y <- seq(-3, 3, len = 1000)
z <- outer(x, y, f)
cr <- contourLines(x, y, z, levels = 0)

library(ggplot2)
dat <- data.frame(x = cr[[1]]$x, y = cr[[1]]$y)
ggplot() + geom_path(aes(x, y), data = dat)

观察到两个现象：

• 直接用contour(x, y, z, levels = 0)能绘制出预期的对称曲线
• 将contourLines中的levels = 0替换为levels = 0.0001可得到对称曲线

请问这是什么原因？如何获取完整的对称曲线？

原因分析 #

1. 等值线追踪的数值敏感性：contourLines的追踪算法对严格等于0的等值面存在精度适配问题。当levels=0时，网格点上恰好等于0的z值分布不均，算法在对称区域的数值变化处容易中断追踪，仅识别出单侧曲线。而contour函数内部有额外的绘图补全逻辑，因此能显示完整对称曲线。
2. 微小偏移的覆盖效应：将levels设为0.0001时，相当于在原曲线的微小邻域内追踪等值线，这个邻域覆盖了对称曲线的两个分支，算法因此能完整识别出全部曲线。

解决方法 #

方法1：利用对称性补全曲线 #

既然曲线关于x轴对称，直接将已获取的单侧曲线的x取反，生成对称分支后合并绘图：

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

x <- seq(-3, 3, len = 1000)
y <- seq(-3, 3, len = 1000)
z <- outer(x, y, f)
cr <- contourLines(x, y, z, levels = 0)

# 提取已获取的单侧曲线
dat <- data.frame(x = cr[[1]]$x, y = cr[[1]]$y)
# 生成对称分支数据
dat_sym <- data.frame(x = -dat$x, y = dat$y)
# 合并完整数据
dat_full <- rbind(dat, dat_sym)

library(ggplot2)
ggplot() + geom_path(aes(x, y), data = dat_full)

方法2：使用接近0的levels值 #

直接设置一个微小的非零levels值，让算法覆盖完整的对称曲线：

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

x <- seq(-3, 3, len = 1000)
y <- seq(-3, 3, len = 1000)
z <- outer(x, y, f)
# 使用接近0的水平值
cr <- contourLines(x, y, z, levels = 1e-4)

library(ggplot2)
# 合并多分支曲线数据
dat_full <- do.call(rbind, lapply(cr, function(line) data.frame(x=line$x, y=line$y)))
ggplot() + geom_path(aes(x, y), data = dat_full)

方法3：使用专门的隐函数绘图工具 #

借助rootSolve包的contourPlot函数，它对隐函数等值线的处理更稳定：

library(rootSolve)
library(ggplot2)

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

# 生成等值线数据但不直接绘图
contour_data <- contourPlot(f, xlim = c(-3, 3), ylim = c(-3, 3), levels = 0, draw = FALSE)
dat_full <- data.frame(x = contour_data$x, y = contour_data$y)

ggplot() + geom_path(aes(x, y), data = dat_full)

内容的提问来源于stack exchange，提问作者Stéphane Laurent