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

Question

我正在尝试找到一种在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, 25,
  'b', 75, 25,
  'c', 50, 75,
) %>%
  {. ->> sample_points}

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

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

expand.grid(
  x = seq(0, 100, length = 100),
  y = seq(0, 100, length = 100)
) %>%
  as_tibble %>%
  mutate(
    # z = f1(x, y, nr)
    z1 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 5),
    z2 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 15),
    z3 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 25),
    z4 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 35),
    z5 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 45),
    z6 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 50),
  ) %>%
  ggplot()+
  geom_contour(aes(x, y, z = z1), breaks = 0, colour = 1)+
  geom_contour(aes(x, y, z = z2), breaks = 0, colour = 2)+
  geom_contour(aes(x, y, z = z3), breaks = 0, colour = 3)+
  geom_contour(aes(x, y, z = z4), breaks = 0, colour = 4)+
  geom_contour(aes(x, y, z = z5), breaks = 0, colour = 5)+
  geom_contour(aes(x, y, z = z6), breaks = 0, colour = 6)+
  geom_point(data = sample_points, aes(x, y, color = name), 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() 之类的函数的选项，但收效甚微，这可能是由于我对基础数学的理解有限和/或描述双曲线的方程同时具有@987654368 @ 和 y 条款？

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

（或三个方程式）：

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

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

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

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

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

Answer 1

更新：

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

查看各种示例：

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

基本上，您为nleqslv() 提供一些起始值* x（来自?nleqslv()：对根的初始猜测）和一个函数fn，其中包含解决。

对于这个例子，我们在问题中指定的z 函数为x 和y 的每个值计算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 等是x 和y 点i 的坐标，r_ij 是交点与点i 和j 之间的距离差.将原方程中的x替换为x[1]，将y替换为x[2]；我的理解是生成的向量（有点混淆，也称为x）将有2个分量，它们是函数fn中的x[1]、x[2]，或原始方程中的x、y（ s) 基于f1()。 x1_start 和 x2_start 是 x[1] (x) 和 x[2] (y) 的起始值。

从我们的问题sample_points，我们分别调用点a, b, c，i, j, k。然后，替换xi, yi, xj, yj, 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}

nleqslv() 生成一个列表，其中包含各种值和消息，解释了方程是如何求解的（有关详细信息，请参阅 ?nleqslv()），但名为 x 的列表的第一个组件包含长度为 2 的向量中的解.

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, 100, length = 100)
) %>%
  as_tibble %>%
  mutate(
    # z = f1(x, y, nr)
    z1 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 5), # i = point a, j = point b
    z2 = f1(xi = 25, yi = 25, xj = 50, yj = 75, x, y, r = 5)  # i = point a, j = point c
  ) %>%
  {. ->> hyp_data} %>%
  ggplot()+
  geom_contour(aes(x, y, z = z1), breaks = 0, colour = 1)+
  geom_contour(aes(x, y, z = z2), breaks = 0, colour = 2)+
  geom_point(data = sample_points, aes(x, y, 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')

所以，我们找到了双曲线的交点——耶！

*选择起始值

但是，起始值的选择仍然存在争议。在这种情况下，我根据sample_points 的最小x 和y 值选择了c(25, 25) 的起始值，因为我“猜到”根（交点）一定在附近。

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

result_list 的一部分是termcd，这是终止代码。据我所见，1 和2 的值表示已找到根，但其他值往往表明未找到（正确的）根（尽管向用户返回了不正确的值）。有关终止代码的完整列表，请参阅?nleqslv（我个人看不懂大部分单词！）。

作为示例，我们使用0 和100 之间的x 和y 的整数组合运行nleqslv()，以查看哪个产生了正确的根，以及返回了termcd 的哪些值：

# set up a grid with every combo of x and y between 0 and 100
expand.grid(
  x_grid = seq(0, 100, length = 101),
  y_grid = seq(0, 100, length = 101)
) %>%
  as_tibble %>%
  rowwise %>%

  # now we run nleqslv(), and extract the x, y, 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, 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]][2],

    termcd = 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
      })[[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, y_grid, label = termcd), size = 2)+
  geom_point(data = sample_points, aes(x, y, color = name), size = 10)+
  geom_text(data = sample_points, aes(x, y, label = name))+
  geom_contour(data = hyp_data, aes(x, y, z = z1), breaks = 0, colour = 1, size = 1)+
  geom_contour(data = hyp_data, aes(x, y, z = z2), breaks = 0, colour = 2, size = 1)+
  labs(
    x = 'starting x',
    y = 'starting y',
    fill = 'Correct answer?',
    colour = 'sample_point name'
    )

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

所以，只要你选择的起始值是正确的，你就应该得到正确的答案。看起来正确答案的终止代码为1 和2，而1 和2 的终止代码只有正确的答案，至少在这种情况下：

# which codes were returned for correct and incorrect answers?
grid_results_2 %>% 
  group_by(termcd) %>% 
  distinct(correct) %>% 
  arrange(termcd)

# # A tibble: 6 x 2
# # Groups:   termcd [6]
#  termcd correct
#   <int> <chr>
# 1     1 Y
# 2     2 Y
# 3     3 N
# 4     4 N
# 5     5 N
# 6     6 N

所以对于这个例子，我猜你可以运行多个选项作为起始值，然后只保留终止代码为1 或2 的值，以确保你得到正确的答案。

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

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