使用SymPy求解Haversine公式坐标分量时程序挂起问题

我碰到过类似的问题，SymPy的符号求解在处理这种带三角函数的超越方程时确实容易卡壳，原因和解决方法我给你梳理一下：

你尝试求解的是一个超越方程（包含三角函数与未知变量的混合表达式），SymPy的solve()函数主打寻找解析解，但这类方程通常没有闭合形式的解析解，SymPy会陷入复杂的符号推导循环，最终导致超时或无限挂起。

因为你要计算的是正北/正南/正东/正西四个特殊方向的点，完全可以简化Haversine公式，直接用数值计算得到结果，比符号求解高效得多：

1. 正北/正南方向（经度不变） #

当两点在同一条经线上时，Haversine公式可以简化为：

球面距离 = 地球半径 × 纬度差（弧度）

推导后，纬度差的弧度值为 Δlat_rad = 距离 / 地球半径，转换为度数后直接加减到中心纬度即可。

2. 正东/正西方向（纬度不变） #

当两点在同一条纬线上时，Haversine公式简化后可以解出经度差：

sin(Δlon/2) = sin(距离/(2×地球半径)) / cos(中心纬度弧度)

计算出经度差的弧度值后，加减到中心经度即可。

代码实现（用Python标准库math） #

import math

# 定义参数
EARTH_RADIUS_MILES = 3950.0
center_lat = 38.0
center_lon = -77.0
target_distance = 5.0

# 计算正北点
delta_lat_rad = target_distance / EARTH_RADIUS_MILES
north_lat = center_lat + math.degrees(delta_lat_rad)
north_point = (round(north_lat, 6), center_lon)

# 计算正南点
south_lat = center_lat - math.degrees(delta_lat_rad)
south_point = (round(south_lat, 6), center_lon)

# 计算正东点
center_lat_rad = math.radians(center_lat)
sin_half_dlon = math.sin(target_distance / (2 * EARTH_RADIUS_MILES)) / math.cos(center_lat_rad)
delta_lon_rad = 2 * math.asin(sin_half_dlon)
east_lon = center_lon + math.degrees(delta_lon_rad)
east_point = (center_lat, round(east_lon, 6))

# 计算正西点
west_lon = center_lon - math.degrees(delta_lon_rad)
west_point = (center_lat, round(west_lon, 6))

# 输出结果
print(f"正北点: {north_point}")
print(f"正南点: {south_point}")
print(f"正东点: {east_point}")
print(f"正西点: {west_point}")

如果你坚持要用SymPy，可以放弃solve()（找解析解），改用nsolve()进行数值求解，它会通过迭代快速找到近似解，不会挂起：

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
R = 3950.0

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))
eq = R * c - d

# 用nsolve，传入初始猜测值（比如中心纬度38.0）
lat2_sol_deg = s.deg(s.nsolve(eq, lat2, 38.0))
print(f"正北点纬度: {round(lat2_sol_deg, 6)}")

这个方法同样适用于其他方向的求解，只需要调整符号变量（比如求解正东点时固定纬度，设经度为变量）和初始猜测值即可。

内容的提问来源于stack exchange，提问作者Brad Solomon