基于直线方程的 R 中双曲线的交点

问题描述

我试图找到一种在 R 中找到 2 个双曲线的交点的方法。单分支双曲线可以用以下等式来描述：

或

其中(xi,yi)和(xj,yj)是2个焦点（i和j）的坐标，r是给定点之间的距离差在双曲线 (x,y) 和每个焦点上。

似乎使用 R 可视化双曲线的最佳方法是可视化 3D 双曲面的轮廓线（当轮廓级别 = 0 时），使用上述等式确定并将其实现为函数

f1 <- function(xi,yi,xj,yj,x,y,r){
  sqrt((xj - x)^2 + (yj - y)^2) - sqrt((xi - x)^2 + (yi - y)^2) - r
  }

例如，我们可以构建一个网格并可视化两个双曲线的 0 级轮廓：

library(tidyverse)
library(sf)

# sample points
tribble(
  ~name,~x,~y,'a',25,'b',75,'c',50,) %>%
  {. ->> sample_points}

# sample hyperbolae
expand.grid(
  x = seq(0,100,length = 100),y = seq(0,length = 100)
) %>%
  as_tibble %>%
  mutate(
    z1 = f1(xi = 25,yi = 25,xj = 75,yj = 25,r = 5),# i = point a,j = point b
    z2 = f1(xi = 25,xj = 50,yj = 75,r = 5)  # i = point a,j = point c
  ) %>%
  {. ->> hyp_data} %>%
  ggplot()+
  geom_contour(aes(x,z = z1),breaks = 0,colour = 1)+
  geom_contour(aes(x,z = z2),colour = 2)+
  geom_point(data = sample_points,aes(x,color = name),size = 3)+
  coord_sf()

目前，我可以通过从两条geom_contour线中提取线数据来找到交点，使用ggplot_build()，将线坐标转换为sf LInesTRING ，并使用 st_intersection 找到交点：

# extract the data from the ggplot using ggplot_build()
ggplot_build(

  hyp_data %>%
    ggplot()+
    geom_contour(aes(x,colour = 1)+
    geom_contour(aes(x,colour = 2)+
    # geom_point(data = sample_points,size = 3)+ # we dont need the point data in ggplot_build()
    coord_sf()

) %>%
  .$data %>%  # keep only the data component
  map(. %>% select(x,y) %>% as.matrix %>% st_linestring) %>% # keep coordinates,turn into a linestring
  {. ->> lines1}

# make the lines an sf object
tibble(a = 1:2) %>%
  mutate(
    geom = st_sfc(lines1),) %>%
  st_as_sf %>% # then make the whole thing an sf object

  # use st_intersection to find the point of intersection
  st_intersection %>%

  # then keep only the 'point' (exclude the original 'lines')
  filter(
    st_is(geom,'POINT') 
  ) %>%
  {. ->> intersection_point} %>%
  ggplot()+
  geom_sf(colour = 'red',size = 5)+
  geom_contour(data = hyp_data,colour = 1)+
  geom_contour(data = hyp_data,colour = name),size = 3)

此方法的局限性在于它依赖于 ggplot 的 geom_contour 可视化轮廓线的能力。当 r 的绝对值接近两点之间的距离（a-b 之间的距离 = 50）时，双曲线变窄并最终成为其中一个点“后面”的直线。此处，geom_contour 无法构建等高线，因此未创建线数据，因此无法找到交点。看看下面的图是如何在应该有 6 条线的情况下只有 5 条线的； z6 不绘图并导致警告消息：

expand.grid(
  x = seq(0,length = 100)
) %>%
  as_tibble %>%
  mutate(
    # z = f1(x,nr)
    z1 = f1(xi = 25,z2 = f1(xi = 25,r = 15),z3 = f1(xi = 25,r = 25),z4 = f1(xi = 25,r = 35),z5 = f1(xi = 25,r = 45),z6 = f1(xi = 25,r = 50),) %>%
  ggplot()+
  geom_contour(aes(x,colour = 2)+
  geom_contour(aes(x,z = z3),colour = 3)+
  geom_contour(aes(x,z = z4),colour = 4)+
  geom_contour(aes(x,z = z5),colour = 5)+
  geom_contour(aes(x,z = z6),colour = 6)+
  geom_point(data = sample_points,size = 3)+
  coord_sf()

# Warning messages:
# 1: stat_contour(): Zero contours were generated 
# 2: In min(x) : no non-missing arguments to min; returning Inf
# 3: In max(x) : no non-missing arguments to max; returning -Inf

z6 理论上应该用下面的虚线表示：

我开发了在无法生成等高线时创建这些直线的方法，但我希望能够使用“数学/代数”方法来找到两条线的交点，而不是依赖 R 中的空间/映射/计数方法。

我已经探索了使用 uniroot() 和 solve() 等函数的选项，但效果有限，可能是因为我对基础数学和/或描述双曲线的方程的理解有限有 x 和 y 项吗？

我目前的想法是写一对等式，右边相等为 0（为不正确的数学语言道歉？）：

（或三个等式）：

其中(xi,yi)、(xj,yj)和(xk,yk)是三个焦点（i、j和k）的坐标，和 r 是范围差异。 xi,xk,yk可以代替上面例子中a,b,c点的坐标

然后，我认为应该可以使用与此处描述的过程类似的过程来求解方程https://www.wikihow.com/Algebraically-Find-the-Intersection-of-Two-Lines，但我不知道这如何转化为 R 代码和/或现有的 R 函数。

这样做的目的是根据接收天线（焦点）接收到的传输的到达时间差 (TDOA) 和双曲线定位原理以及三边测量/多边测量来确定发射机的位置（交点） （类似于 GPS 的工作原理）。

我将处理数千到数百万对双曲线，因此理想情况下，此过程相对较快且可由普通计算机管理。此外，不是必需的，但希望可以实现是在相同传输被三个以上的电台（例如使用 GPS 时的“卫星数量”）。

我非常感谢任何人可以提供的任何想法或帮助，因为我一直在思考并努力解决这个问题已经有一段时间了。谢谢。

解决方法

更新：

对于在家玩耍的任何人来说，nleqslv::nleqslv() 函数都可用于求解联立方程（非线性方程求解器；n-l-eq-slv()）。

查看各种示例：

• https://stackoverflow.com/a/48833164/13478749
• Solving system of nonlinear equations in R

基本上，您为 nleqslv() 提供一些起始值* x（来自 ?nleqslv()：对根的初始猜测）和一个函数 {{ 1}} 包含要求解的方程。

对于这个例子，我们在问题中指定的 fn 函数计算 f1() 和 z 的每个值的 x 值，是每个双曲线的方程，并形成待解方程。

y

当使用 f1 <- function(xi,yi,xj,yj,x,y,r){ sqrt((xj - x)^2 + (yj - y)^2) - sqrt((xi - x)^2 + (yi - y)^2) - r } 时，这会转化为类似的内容：

nleqslv()

其中nleqslv::nleqslv(c(x1_start,x2_start),function(x){ z <- numeric() z[1] = sqrt((xj - x[1])^2 + (yj - x[2])^2) - sqrt((xi - x[1])^2 + (yi - x[2])^2) - r_ij z[2] = sqrt((xk - x[1])^2 + (yk - x[2])^2) - sqrt((xi - x[1])^2 + (yi - x[2])^2) - r_ik z }) 、xi等是点yi的{​​{1}}和x坐标，y是点之间的距离差交点和点 i 和 r_ij。原方程中的i替换为j，x替换为x[1]；我的理解是生成的向量（也命名为 y 有点令人困惑）将有 2 个分量，它们是函数 x[2] 中的 x、x[1] 或 {{1} },x[2] 在基于 fn 的原始方程中。 x 和 y 是 f1() (x1_start) 和 x2_start (x[1]) 的起始值。

从问题的 x 中，我们分别称点 x[2]、y。然后，代入 sample_points 的值和 a,b,c 的一些起始值*，我们得到如下所示的结果：

i,j,k

xi,xk,yk,d_ij,d_ik 生成一个列表，其中包含解释方程如何求解的各种值和消息（有关详细信息，请参阅 c(25,25)），但名为 library(nleqslv) nleqslv::nleqslv(c(25,25),function(x){ z <- numeric() z[1] = sqrt((75 - x[1])^2 + (25 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z[2] = sqrt((50 - x[1])^2 + (75 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z }) %>% {. ->> result_list} 的列表的第一个组件包含解决方案，在长度为 2 的向量中。

nleqslv()

第一个元素是我们对 ?nleqslv()（原方程中的x）的值，第二个元素是 result_list[[1]] # [1] 46.95886 42.22943 （x[1]）。

如我们所见，该点与我们原始双曲线的交点相匹配：

x

所以，我们找到了双曲线的交点 - 是的！

*选择起始值

然而，起始值的选择仍然存在争议。在本例中，我根据 x[2] 的最小 y 和 # extract the values of x and y from result_list for plotting px <- result_list[[1]][1] py <- result_list[[1]][2] # plot the sample hyperbolae as in the question expand.grid( x = seq(0,100,length = 100),y = seq(0,length = 100) ) %>% as_tibble %>% mutate( # z = f1(x,nr) z1 = f1(xi = 25,yi = 25,xj = 75,yj = 25,r = 5),# i = point a,j = point b z2 = f1(xi = 25,xj = 50,yj = 75,r = 5) # i = point a,j = point c ) %>% {. ->> hyp_data} %>% ggplot()+ geom_contour(aes(x,z = z1),breaks = 0,colour = 1)+ geom_contour(aes(x,z = z2),colour = 2)+ geom_point(data = sample_points,aes(x,color = name),size = 3)+ coord_sf()+ # and add in the newfound point geom_point(aes(x = px,y = py),size = 5,shape = 1,colour = 'red') 值选择了 c(25,25) 的起始值，因为我“猜测”了根（交点）必须在这附近的某个地方。

请注意，替代起始值可能会返回不正确的根。幸运的是，似乎有一些方法可以尝试解决这个问题。

x 的一部分是 y，它是终止代码。据我所知，sample_points 和 result_list 的值表示已找到根，但其他值往往表示未找到（正确的）根（尽管向用户返回了不正确的值） .请参阅 termcd 以获取终止代码的完整列表（我个人不理解大部分单词！）。

例如，我们使用 1 和 2 之间的 ?nleqslv 和 nleqslv() 的整数组合运行 x 以查看哪个产生正确的根，并且返回 y 的哪些值：

0

蓝色单元格是返回正确根的起始值，红色单元格返回错误的根。数字是 100 值，蓝色区域是 1 或 2，红色区域是其他值。

所以，只要你选择的起始值是正确的，你就应该得到正确的答案。看起来正确的答案只有 termcd 和 # set up a grid with every combo of x and y between 0 and 100 expand.grid( x_grid = seq(0,length = 101),y_grid = seq(0,length = 101) ) %>% as_tibble %>% rowwise %>% # now we run nleqslv(),and extract the x,and termcd values mutate( px = nleqslv::nleqslv(c(x_grid,y_grid),function(x){ z <- numeric() z[1] = sqrt((75 - x[1])^2 + (25 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z[2] = sqrt((50 - x[1])^2 + (75 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z })[[1]][1],py = nleqslv::nleqslv(c(x_grid,function(x){ z <- numeric() z[1] = sqrt((75 - x[1])^2 + (25 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z[2] = sqrt((50 - x[1])^2 + (75 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z })[[1]][2],termcd = nleqslv::nleqslv(c(x_grid,function(x){ z <- numeric() z[1] = sqrt((75 - x[1])^2 + (25 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z[2] = sqrt((50 - x[1])^2 + (75 - x[2])^2) - sqrt((25 - x[1])^2 + (25 - x[2])^2) - 5 z })[[3]][1],) %>% {. ->> grid_results} # in this instance,we know the root/intersection is roughly (47,42),so let's check # which starting values produce the correct result grid_results %>% ungroup %>% # work out which starting values got the correct point mutate( correct = ifelse(round(px) == 47 & round(py) == 42,'Y','N') ) %>% {. ->> grid_results_2} %>% # then plot,showing colour for correct/incorrect,and include termcd as text ggplot()+ geom_raster(aes(x_grid,y_grid,fill = factor(correct)))+ coord_sf()+ geom_text(aes(x_grid,label = termcd),size = 2)+ geom_point(data = sample_points,size = 10)+ geom_text(data = sample_points,label = name))+ geom_contour(data = hyp_data,colour = 1,size = 1)+ geom_contour(data = hyp_data,colour = 2,size = 1)+ labs( x = 'starting x',y = 'starting y',fill = 'Correct answer?',colour = 'sample_point name' ) 的终止代码，而 termcd 和 1 的终止代码只有正确的答案，至少在这种情况下是这样：

2

因此，对于此示例，我猜您可以针对起始值运行多个选项，然后仅保留终止代码为 1 或 2 的值，以确保您得到正确的答案。

我想这种方法的一个限制是您必须为每个实例运行 10,000 次初始值迭代 - 当有 10 到 100 个数千（或数百万）个实例时，这可能会减慢处理时间。为了缓解这种情况，我使用 25 个起始值而不是 10,000 运行了与上面相同的代码，并且在大多数情况下它返回了正确的答案，终止代码的模式相同：

目前，这是一个可行的解决方案。