如何用Python求解含X、Y的二元三角函数方程组

这组方程看起来是典型的球面坐标转换相关问题，我会先帮你推导解析解，再给出Python的实现代码，同时也会介绍数值解法作为备选方案——毕竟有时候手动推导容易漏特殊情况，数值法反而更省心。

一、数学推导（解析解） #

先做个变量替换简化问题：令 $\theta = Y - b$，这样原方程组可以写成：

1. $\cos(z) = \sin(a)\sin(X) + \cos(a)\cos(X)\cos(\theta)$
2. $\tan(Az) = \frac{\cos(X)\sin(\theta)}{\cos(a)\sin(X) - \sin(a)\cos(X)\cos(\theta)}$

通过三角函数恒等变换（核心是利用$\sin^2x+\cos2x=1$和$\tan^2x+1=\sec2x$消元），最终可以得到两组配对的解析解：

解的表达式 #

两组解的符号要保持对应：

1. 计算$X$：
   $\sin(X) = \sin(a)\cos(z) \pm \cos(a)\sin(z)\cos(Az)$
   $X = \arcsin\left(\sin(a)\cos(z) \pm \cos(a)\sin(z)\cos(Az)\right)$
   （注意$\arcsin$返回范围是$[-\pi/2, \pi/2]$，如果你的实际问题中$X$有其他范围要求，需要额外调整）

2. 计算$\theta$：
   $\theta = \arctan2\left(\pm \sin(z)\sin(Az), \cos(a)\cos(z) \mp \sin(a)\sin(z)\cos(Az)\right)$

3. 计算$Y$：
   $Y = \theta + b$

这里的$\pm$和$\mp$是配对的：取上面的符号时，所有$\pm$用正号、$\mp$用负号；取下面的符号时则相反，对应两组不同的解。

二、Python解析法实现 #

注意：所有角度参数必须用弧度制，如果你的输入是角度，记得先用math.radians()转换。

import math

def solve_spherical_equations(z_rad, Az_rad, a_rad, b_rad):
    solutions = []
    # 遍历两组符号组合
    for sign in [1, -1]:
        # 计算sin(X)，并限制在[-1,1]范围内（避免浮点误差导致超出）
        sin_X = math.sin(a_rad) * math.cos(z_rad) + sign * math.cos(a_rad) * math.sin(z_rad) * math.cos(Az_rad)
        sin_X = max(min(sin_X, 1.0), -1.0)
        X_rad = math.asin(sin_X)
        
        # 计算theta的分子项
        cos_theta_num = math.cos(a_rad) * math.cos(z_rad) - sign * math.sin(a_rad) * math.sin(z_rad) * math.cos(Az_rad)
        sin_theta_num = sign * math.sin(z_rad) * math.sin(Az_rad)
        
        # 计算theta和Y
        theta_rad = math.atan2(sin_theta_num, cos_theta_num)
        Y_rad = theta_rad + b_rad
        
        # 验证解是否满足原方程（可选，排查数值误差）
        def verify(X, Y):
            calc_z = math.acos(math.sin(a_rad)*math.sin(X) + math.cos(a_rad)*math.cos(X)*math.cos(Y - b_rad))
            calc_Az = math.atan2(math.cos(X)*math.sin(Y - b_rad), math.cos(a_rad)*math.sin(X) - math.sin(a_rad)*math.cos(X)*math.cos(Y - b_rad))
            return math.isclose(calc_z, z_rad, abs_tol=1e-6) and math.isclose(calc_Az, Az_rad, abs_tol=1e-6)
        
        if verify(X_rad, Y_rad):
            solutions.append((X_rad, Y_rad))
    
    return solutions

# 示例用法
if __name__ == "__main__":
    # 假设已知参数（先转成弧度）
    z = math.radians(30)
    Az = math.radians(45)
    a = math.radians(60)
    b = math.radians(15)
    
    results = solve_spherical_equations(z, Az, a, b)
    print("所有有效解（转换为角度）：")
    for i, (X, Y) in enumerate(results):
        print(f"解{i+1}: X={math.degrees(X):.2f}°, Y={math.degrees(Y):.2f}°")

三、数值解法（无需手动推导） #

如果解析推导太繁琐，或者担心特殊情况（比如$\cos(a)=0$或$\cos(X)=0$）处理不到位，可以用数值优化的方法，直接调用scipy.optimize.root求解残差函数：

from scipy.optimize import root
import math

def residual(vars, z_rad, Az_rad, a_rad, b_rad):
    X, Y = vars
    # 计算两个方程的残差（即左边减右边）
    res1 = math.acos(math.sin(a_rad)*math.sin(X) + math.cos(a_rad)*math.cos(X)*math.cos(Y - b_rad)) - z_rad
    res2 = math.atan2(math.cos(X)*math.sin(Y - b_rad), math.cos(a_rad)*math.sin(X) - math.sin(a_rad)*math.cos(X)*math.cos(Y - b_rad)) - Az_rad
    return [res1, res2]

def solve_numerically(z_rad, Az_rad, a_rad, b_rad, initial_guess=[0, 0]):
    result = root(residual, initial_guess, args=(z_rad, Az_rad, a_rad, b_rad))
    if result.success:
        X, Y = result.x
        return [(X, Y)]
    else:
        print("数值求解失败：", result.message)
        return []

# 示例用法
if __name__ == "__main__":
    z = math.radians(30)
    Az = math.radians(45)
    a = math.radians(60)
    b = math.radians(15)
    
    num_results = solve_numerically(z, Az, a, b)
    print("数值解（转换为角度）：")
    for X, Y in num_results:
        print(f"X={math.degrees(X):.2f}°, Y={math.degrees(Y):.2f}°")

注意事项 #

• 解析法会得到两组解，需要根据你的实际业务场景（比如物理意义、角度范围）筛选符合要求的解；
• 数值解法依赖初始猜测，如果初始值离真实解太远可能收敛失败，建议根据问题范围设置合理的初始值；
• 所有三角函数计算都使用弧度制，角度输入务必转换；
• 浮点计算存在微小误差，验证时用math.isclose而非直接相等判断。

内容的提问来源于stack exchange，提问作者Neeraj Sirdeshmukh