如何使用 SymPy 求解方程组来找到一些角度？

Question

我不太明白如何在 Python 中使用 SymPy 求解方程组。我尝试过使用 solve() 和 nsolve()。 solve() 没有给我答案，它只是静止不动，nsolve() 给了我这个错误。

ValueError：无法在给定的容差范围内找到根。 (7378.8100001153525709 > 2.16840434497100886801e-19) 尝试另一个起点或调整参数。

这是我在尝试使用 solve() 时所写的内容：

from sympy import *

theta_1, theta_2, theta_3 = symbols('theta_1 theta_2 theta_3')

eq1 = -136.2*sin(theta_2)*sin(theta_3)*cos(theta_1) + 136.2*cos(theta_1)*cos(theta_2)*cos(theta_3) + 222.1*cos(theta_1)*cos(theta_2)

eq2 = -136.2*sin(theta_1)*sin(theta_2)*cos(theta_3) + 136.2*sin(theta_1)*cos(theta_2)*cos(theta_3) + 222.1*sin(theta_1)*cos(theta_2)

eq3 =  136.2*sin(theta_2)*cos(theta_3)+ 222.1*sin(theta_2) + 136.2*sin(theta_3)*cos(theta_2) + 100.9 

eq1 = Eq(eq1, 0)
eq2 = Eq(eq2, (-323.90))
eq3 = Eq(eq3, 176.71)

ans = solve((eq1, eq2, eq3), (theta_1, theta_2, theta_3))

print(ans)

这是我在尝试使用 solve() 时所写的内容 ：尝试使用nsolve()：

from sympy import *

theta_1, theta_2, theta_3 = symbols('theta_1 theta_2 theta_3')

eq1 = -136.2*sin(theta_2)*sin(theta_3)*cos(theta_1) + 136.2*cos(theta_1)*cos(theta_2)*cos(theta_3) + 222.1*cos(theta_1)*cos(theta_2)

eq2 = -136.2*sin(theta_1)*sin(theta_2)*cos(theta_3) + 136.2*sin(theta_1)*cos(theta_2)*cos(theta_3) + 222.1*sin(theta_1)*cos(theta_2)

eq3 =  136.2*sin(theta_2)*cos(theta_3)+ 222.1*sin(theta_2) + 136.2*sin(theta_3)*cos(theta_2) + 100.9 

ans = nsolve((eq1, eq2, eq3), (theta_1, theta_2, theta_3), (0, -323.9, 176.71))

print(ans)

我想找到一个表达式θ1、θ2 和 θ3。优选为 arctan2 的形式。有谁知道我该怎么做？或者为什么这不起作用？

Answer 1

我假设您的输入应该是以度为单位的角度。如果您以弧度为单位给出输入，nsolve 可以给您答案：

In [48]: p = (pi/180).evalf()

In [49]: nsolve((eq1, eq2, eq3), (theta_1, theta_2, theta_3), (0, -323.9*p, 176.71*p))
Out[49]: 
⎡-1.5707963267949 ⎤
⎢                 ⎥
⎢-9.36640725122338⎥
⎢                 ⎥
⎣-2.49008907692194⎦

我们可以将您的系统简化为具有有理系数的多项式系统，如下所示：

In [60]: c1, c2, c3, s1, s2, s3 = symbols('c1:4 s1:4')
    ...: rep = {cos(theta_1): c1, cos(theta_2): c2, cos(theta_3): c3,
    ...:        sin(theta_1): s1, sin(theta_2): s2, sin(theta_3): s3}
    ...: 
    ...: eqs = [eq.subs(rep) for eq in [eq1, eq2, eq3]] + [
    ...:         c1**2 + s1**2 - 1, c2**2 + s2**2 - 1, c3**2 + s3**2 - 1]
    ...: eqs = [nsimplify(eq) for eq in eqs]

In [61]: for eq in eqs:
    ...:     print(eq)
    ...: 
Eq(681*c1*c2*c3/5 - 681*c1*s2*s3/5 + 2221*c1/10, 0)
Eq(681*c2*c3*s1/5 - 681*c3*s1*s2/5 + 2221*s1/10, -3239/10)
Eq(681*c2*s3/5 + 681*c3*s2/5 + 2221*s2/10 + 1009/10, 17671/100)
c1**2 + s1**2 - 1
c2**2 + s2**2 - 1
c3**2 + s3**2 - 1

从那里您可以使用 SymPy 的 groebner 函数计算多项式系统解的 Groebner 基。为了对称性，我们将添加一个分隔元素 t：

In [62]: gb = groebner(eqs, order='grevlex').fglm('lex')

这就像一个新的分隔方程列表，其解是您想要的不同 theta 的 cos 和 sin 值。

最终方程是 c3 的 20 次多项式，其解可以代入其他方程以给出系统的完整解。该多项式有 4 个实根，每个实根都是一个较小的 8 次不可约多项式的根。由于 Abel-Ruffini，根不能有根式表达式，但可以用数值方法计算它们：

In [70]: c3num = [r.n() for r in real_roots(gb[-1])]

In [71]: c3num
Out[71]: [-0.795172954862362, -0.552080421425289, 0.702791508766261, 0.99923774387104]

其中每一个都可用于求解 Groebner 基的其余部分：

In [73]: for c3n in c3num:
     ...:     for eq in gb[:-1]:
     ...:         print(solve([eq.subs(c3, c3n)], dict=True))
     ...:     print(solve([c3 - c3n], dict=True))
     ...:     print()
     ...: 
[{s1: -1.00000000000000}]
[{s2: -0.0583375690202956}]
[{s3: -0.606382694211788}]
[{c1: 0}]
[{c1: 0.0}]
[{c2: -0.998296913792728}]
[{c3: -0.795172954862362}]

[{s1: -1.00000000000000}]
[{s2: 0.881316344165583}]
[{s3: 0.833790866034178}]
[{c1: 0}]
[{c1: 0.0}]
[{c2: -0.472526720395763}]
[{c3: -0.552080421425289}]

[{s1: -1.00000000000000}]
[{s2: -0.0656752871541357}]
[{s3: 0.711395878031908}]
[{c1: 0}]
[{c1: 0.0}]
[{c2: 0.997841047763359}]
[{c3: 0.702791508766261}]

[{s1: -1.00000000000000}]
[{s2: 0.226102985521720}]
[{s3: -0.0390375619754195}]
[{c1: -3.05175781250000e-5}, {c1: 3.05175781250000e-5}]
[{c1: 0.0}]
[{c2: 0.974103403161280}]
[{c3: 0.999237743871040}]

最后一种情况似乎在数值上不一致，但这可能是舍入误差，因为这在数值上相当不稳定，也许解应该是 c1=0。我认为发生这种情况是因为你的方程是一个稍微退化的情况。

上面给出了角度的余弦和正弦的所有解决方案，您可以使用 atan2 从那里得到角度本身。您可以使用 nsolve 根据需要精确地优化上述任何内容。 4 个案例中的第一个是上面 nsolve 发现的案例。

如果你想要一个明确的分析表达式，那么这似乎不太可能。