sympy：如何解决这个超越方程组？

我正试图用sympy来解决一个超越方程组：

from sympy import *
A, Z, phi, d, a, tau, t, R, u, v = symbols('A Z phi d a tau t R u v', real=True)

system = [
    Eq(A, sqrt(R**2*cos(u)**2 + 2*R*d*cos(a)*cos(u) + d**2 + 2*d*(v*sin(tau) - (R*sin(u) + t)*cos(tau))*sin(a) + (v*sin(tau) - (R*sin(u) + t)*cos(tau))**2)),
    Eq(Z, v*cos(tau) + (R*sin(u) + t)*sin(tau)),
    Eq(phi, (R*sin(a)*cos(u) - (v*sin(tau) - (R*sin(u) + t)*cos(tau))*cos(a))/A),
]

或在欧几里德坐标中：

system = [
     Eq(A, sqrt(d**2 + 2*d*x*cos(a) + 2*d*(z*sin(tau) - (t + y)*cos(tau))*sin(a) + x**2 + (z*sin(tau) - (t + y)*cos(tau))**2)),
     Eq(Z, z*cos(tau) + (t + y)*sin(tau)),
     Eq(phi, (x*sin(a) - (z*sin(tau) - (t + y)*cos(tau))*cos(a))/A),
     Eq(R**2, x**2+y**2),
]

基本上，这是圆（半径A）与围绕偏心轴（d）旋转（角度α）的倾斜（角τ）圆柱（半径R）的交点 . 我对函数φ（A，Z）感兴趣，因为我需要对其他参数（R，A，d，t，τ）的几种组合的A和Z的几百种组合进行评估 .

这是我尝试过的：

• 直接使用 solve(system, [u, v, phi]) ：需要太多时间才能完成 .

• 先使用 solveset() 将其缩小为一个等式 .

• 为所有参数插入数字（例如： R, A, d, t, tao = 25, 150, 125, 2, 5/180*pi ），然后使用 N() 和 solve() 或 solveset() 的任意组合 .

• 我能够使用 d = sqrt(A**2 - R**2) 获得解决了未被歪曲的圆柱体（ t, tao = 0, 0 ）的简化问题 .

How can I solve this system of equations? Failing that: How can I get numeric values for φ for a few dozen combinations of the parameters and a few hundred points (A, Z) each?

如果你感兴趣 . 以下是导致等式的代码：

A, Z, phi, d, a, tau, t, R, u, v = symbols('A Z phi d a tau t R u v', real=True)

r10, z10, phi10 = symbols('r10 z10 phi_10', real=True)
system10 = [Eq(A, r10), Eq(Z, z10), Eq(phi, phi10)]

x9, x8, y8, z8 = symbols('x9 x8 y8 z8', real=True)
system8 = [e.subs(r10, sqrt(x9**2+y8**2)).subs(z10, z8).subs(phi10, y8/A).subs(x9, x8+d) for e in system10]

r6, phi7, phi6, z6 = symbols('r6 phi_7 phi_6 z6', real=True)
system6 = [e.subs(x8, r6*cos(phi7)).subs(y8, r6*sin(phi7)).subs(z8, z6).subs(phi7, phi6+a) for e in system8]

x6, y6 = symbols('x_6 y_6', real=True)
system6xy = [simplify(expand_trig(e).subs(cos(phi6), x6/r6).subs(sin(phi6), y6/r6)) for e in system6]

# solve([simplify(e.subs(phi6, u).subs(r6, R).subs(z6, v).subs(d, sqrt(A**2-R**2))) for e in system6], [u, v, phi]) -> phi = acos(+-R/A)

r5, phi5, x5 = symbols('r_5 phi_5 x_5', real=True)
system5 = [e.subs(x6, x5).subs(y6, r5*cos(phi5)).subs(z6, r5*sin(phi5)) for e in system6xy]

r4, phi4, x4 = symbols('r_4 phi_4 x_4', real=True)
system4 = [e.subs(phi5, phi4 + tau).subs(r5, r4).subs(x5, x4) for e in system5]

x3, y3, z3 = symbols('x_3 y_3 z_3', real=True)
system3 = [expand_trig(e).subs(cos(phi4), y3/r4).subs(sin(phi4), z3/r4).subs(r4, sqrt(y3**2+z3**2)).subs(x4, x3) for e in system4]

x2, y2, z2 = symbols('x_2 y_2 z_2', real=True)
system2 = [e.subs(x3, x2).subs(y3, y2+t).subs(z3, z2) for e in system3]

system = [simplify(e.subs(x2, R*cos(u)).subs(y2, R*sin(u)).subs(z2, v)) for e in system2]

# solve([simplify(e.subs(tau, 0).subs(t, 0).subs(d, sqrt(A**2-R**2))) for e in system], [u, v, phi]) -> phi = acos(+-R/A)