求解Haversine公式的分量时，sympy挂起[重复]

Posted 2023-02-23

【中文标题】求解Haversine公式的分量时，sympy挂起[重复]

【英文标题】：sympy hangs when solving for component of Haversine formula [duplicate]

【发布时间】：2020-07-18 22:23:50

【问题描述】：

我正在使用Haversine formula 计算两个(lat, lng) 坐标对之间的距离（以英里为单位）。 （注意：我知道并且可以接受与地球的非球形（椭圆体）形状相关的公式的限制。）

我想使用这个公式来求解一个纬度或一个经度，它是给定坐标的正北、正东、正西或正南。这可能最好通过图表来说明； 我有给定的中心红点，并试图解决下面的 4 个外部红点：

从 (38.0, -77.0) 的中心坐标开始，假设距离为 5 英里，我想 (单独) 求解图中圆每一侧的 4 个缺失点。所以在每个方程中，给定一个距离和 3 个坐标，并想求解第 4 个坐标。

我尝试使用sympy，但计算似乎超时，除非我在某处有错误的符号。以顶点（lat2, -77.0)为例：

import sympy as s

lat1 = s.rad(38.0)
lat2 = s.Symbol('lat2')
lon1 = s.rad(-77.0)
lon2 = s.rad(-77.0)
d = 5.0  # Given distance
R = 3950.  # Radius of earth in miles
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (s.sin(dlat/2))**2 + s.cos(lat1) * s.cos(lat2) * (s.sin(dlon/2))**2
c = 2 * s.atan2( s.sqrt(a), s.sqrt(1-a) )
s.solve(3950 * c - d, lat2)  # HANGS

我正在尝试解决lat2，但.solve() 呼叫无限期挂断。

【问题讨论】：

Python 中的幂运算符是**，而^ 是异或。

对不起，我知道@MichaelButscher ;) 这里复制了错误的公式。我已经更正了，但我看到了同样的问题。

这能回答你的问题吗？ ***.com/questions/15258078/…

很可能@Joni，但在我们关闭为重复项之前，让我更彻底地了解一下

【参考方案1】：

由于方程中有浮点数并且没有符号，因此您无法获得精确的解析解，因此您可以使用nsolve。请注意，对于不同的初始估计，我看到了不同的解决方案，它们似乎都是解决方案。我不确定如何根据您的问题来解释这一点（也许是 4 点中的 2 点？）。

In [25]: nsolve(eq, lat2, 0)                                                                                           
Out[25]: 0.661959292973035

In [26]: nsolve(eq, lat2, 1)                                                                                           
Out[26]: 0.664490938542655

In [27]: v1 = nsolve(eq, lat2, 0)                                                                                      

In [28]: v2 = nsolve(eq, lat2, 2)                                                                                      

In [29]: eq.subs(lat2, v1).n()                                                                                         
Out[29]: -7.37924650738106e-14

In [30]: eq.subs(lat2, v2).n()                                                                                         
Out[30]: 1.27083170255818e-13

如果你使用 atan 而不是 arctan2：

c = 2 * s.atan( s.sqrt(a) / s.sqrt(1-a) ) 

那么你可以直接从solve获得这两个：

In [43]: s.solve(3950 * c - d, lat2, check=False, simplify=False)                                                      
Out[43]: [0.661959292973035, 0.664490938542655]

我不确定 atan2 出了什么问题...