Как найти пересечение кривой и круга?

У меня есть кривая, полученная из эмпирических данных, и я могу получить разумную модель. Мне нужно определить точку (x, y), где кривая пересекает окружность известного центра и радиуса. Следующий код иллюстрирует вопрос.

x <- c(0.05, 0.20, 0.35, 0.50, 0.65, 0.80, 0.95, 
   1.10, 1.25, 1.40, 1.55, 1.70, 1.85, 2.00, 
   2.15, 2.30, 2.45, 2.60, 2.75, 2.90, 3.05)

y <- c(1.52, 1.44, 1.38, 1.31, 1.23, 1.15, 1.06,
   0.96, 0.86, 0.76, 0.68, 0.61, 0.54, 0.47, 
   0.41, 0.36, 0.32, 0.29, 0.27, 0.26, 0.26)

fit <- loess(y ~ x, control = loess.control(surface = "direct"))
newx <- data.frame(x = seq(0, 3, 0.01))
fitline <- predict(fit, newdata = newx)
est <- data.frame(newx, fitline)

plot(x, y, type = "o",lwd = 2)
lines(est, col = "blue", lwd = 2)

library(plotrix)
draw.circle(x = 3, y = 0, radius = 2, nv = 1000, lty = 1, lwd = 1)

enter image description here

Ответы

Ответ 1

Прямо найти пересечение, используя функции из пакета sf.

Вычислите значения круга (вдохновленные этим ответом и сделанные @Onyambu)

circ <- function(xc = 0, yc = 0, r = 1, n = 100){
  v <- seq(0, 2 * pi, len = n)
  cbind(x = xc + r * cos(v),
        y = yc + r * sin(v))
}

m <- circ(xc = 3, yc = 0, r = 2)

Преобразуйте предсказанные значения и значения круга в "простые функции" (LINESTRING) и найдите их пересечение (POINT):

library(sf)
int <- st_intersection(st_linestring(as.matrix(est)),
                       st_linestring(m))
int
# POINT (1.2091 0.8886608)

Добавьте пересечение к вашему сюжету:

plot(x, y, type = "o", lwd = 2)
lines(est, col = "blue", lwd = 2)
lines(m)
points(int[1], int[2], col = "red", pch = 19)

enter image description here

Ответ 2

Чтобы получить точку пересечения, мы можем использовать функцию оптимизации в г:

circle=function(x){
  if(4<(x-3)^2) return(NA)# Ensure it is limited within the radius
  sqrt(4-(x-3)^2)
}
fun=function(x)predict(fit,data.frame(x=x))  
g=function(x)(circle(x)-fun(x))# We need to set this to zero. Ie solve this
sol1=optimise(function(x)abs(g(x)),c(1,5))$min
 [1] 1.208466

Таким образом, две функции должны оцениваться с одинаковым значением при x=1.208466..

Чтобы сделать его более точным, вы можете использовать функцию optim:

sol2= optim(1,function(x)abs(g(x)),g,method="Brent",upper=5,lower=1)$par
 [1] 1.208473

Теперь вы можете оценить:

circle(sol1)
[1] 0.889047
fun(sol1)
        1 
0.8890654 
circle(sol2)
[1] 0.889061
fun(sol2)
       1 
0.889061

Из вышесказанного вы можете сказать, что решение 2 очень близко.

Построение этой точки на графике будет сложным, так как функция draw.circle рисует круги пропорционально с zxes.. Таким образом, меняется каждый раз в зависимости от того, насколько велика область сюжета.

Если вы должны написать свою собственную функцию круга:

circleplot=function(x,y,r){
  theta=seq(0,2*pi,length.out = 150)
  cbind(x+r*cos(theta),y+r*sin(theta))
}

Тогда вы можете сделать:

plot(x, y, type = "o",lwd = 2)
lines(est, col = "blue", lwd = 2)
lines(circleplot(3,0,2))
abline(v=sol2,col=2) 
points(sol2,fun(sol2),col=2,pch=16)

enter image description here