如何解决基于直线方程的 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
我开发了在无法生成等高线时创建这些直线的方法,但我希望能够使用“数学/代数”方法来找到两条线的交点,而不是依赖 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()
)。
查看各种示例:
基本上,您为 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 运行了与上面相同的代码,并且在大多数情况下它返回了正确的答案,终止代码的模式相同:
目前,这是一个可行的解决方案。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。