python - sympy:双摆方程组

Question

我想求解文章 Double Pendulum 中描述的方程组，与作者类似，在最后一步中他声称他“使用计算机代数程序来求解方程(13)和(16)”。为此，我在 SymPy 中重写了这个方程:

import sympy as s
m1, m2 = s.symbols('m1, m2')
g = s.symbols('g')
x1_d2, y1_d2, x2_d2, y2_d2 = s.symbols('x1_d2, y1_d2, x2_d2, y2_d2')
t1, t2, t1_d1, t1_d2, t2_d1, t2_d2 = s.symbols('t1, t2, t1_d1, t1_d2, t2_d1, t2_d2')
L1, L2 = s.symbols('L1, L2')

eq1 = x1_d2 + t1_d1**2 * L1 * s.sin(t1) - t1_d2 * L1 * s.cos(t1)
eq2 = y1_d2 - t1_d1**2 * L1 * s.cos(t1) - t1_d2 * L1 * s.sin(t1)
eq3 = x2_d2 - x1_d2 + t2_d1**2 * L2 * s.sin(t2) - t2_d2 * L2 * s.cos(t2)
eq4 = y2_d2 - y1_d2 - t2_d1**2 * L2 * s.cos(t2) - t2_d2 * L2 * s.sin(t2)

eq16 = s.sin(t2) * (m2 * y2_d2 + m2 * g) + s.cos(t2) * (m2 * x2_d2)
result1 = s.solve([eq1, eq2, eq3, eq4, eq16], \
    [m1, m2, g, x1_d2, y1_d2, x2_d2, y2_d2, t1, t2, t1_d1, t1_d2, t2_d1, t2_d2, L1, L2], \
    dict=True)

for r in result1:
    print(r.keys())

eq13 = s.sin(t1) * (m1 * y1_d2 + m2 * y2_d2 + m2 * g + m1 * g) + s.cos(t1) * (m1 * x1_d2 + m2 * x2_d2)
result2 = s.solve([eq1, eq2, eq3, eq4, eq13], \
    [m1, m2, g, x1_d2, y1_d2, x2_d2, y2_d2, t1, t2, t1_d1, t1_d2, t2_d1, t2_d2, L1, L2], \
    dict=True)

for r in result2:
    print(r.keys())

正在寻找类似的问题，如 How can I solve system of linear equations in SymPy? ，我期望 SymPy 简化并返回每个符号的方程，但是对于方程 (16)，我仅得到以下结果:L1, g, t1, L2, t2，而不是 y2_d2 t2_d2。对于方程(13)我有异常(exception)。

dict_keys([L1, g, t1, L2, t2])
dict_keys([L1, g, t1, L2, t2])

Traceback (most recent call last):
  File "physics-simulations/formula.py", line 28, in <module>
    dict=True)
  File "/usr/lib/python3/dist-packages/sympy/solvers/solvers.py", line 1164, in solve
    solution = _solve_system(f, symbols, **flags)
  File "/usr/lib/python3/dist-packages/sympy/solvers/solvers.py", line 1911, in _solve_system
    soln = _solve(eq2, s, **flags)
  File "/usr/lib/python3/dist-packages/sympy/solvers/solvers.py", line 1752, in _solve
    result = [r for r in result if
  File "/usr/lib/python3/dist-packages/sympy/solvers/solvers.py", line 1753, in <listcomp>
    checksol(f_num, {symbol: r}, **flags) is not False]
  File "/usr/lib/python3/dist-packages/sympy/solvers/solvers.py", line 355, in checksol
    return bool(abs(val.n(18).n(12, chop=True)) < 1e-9)
  File "/usr/lib/python3/dist-packages/sympy/core/expr.py", line 336, in __lt__
    raise TypeError("Invalid NaN comparison")
TypeError: Invalid NaN comparison

编辑: 我的代码应该是什么样子来求解/找到两个未知数 y1_d2, y2_d2 t1_d1, t2_d2 - 方程 (13, 16)？

Answer 1

描述了这种操作 here ;使用此处描述的 focus 例程，您可以将完整的方程组(转换为 Eq 实例)提供给 focus 并指定 y1_d2 和 y2_d2 作为您要关注的变量:

>>> F = (y1_d2, y2_d2)
>>> f = focus([Eq(i,0) for i in [eq1, eq2, eq3, eq4, eq16, eq13]], *F)
>>> [f[i].has(*F) for i in F]
[False, False]
>>> count_ops(f)
323

使用cse可以获得紧凑的结果:

>>> r, e =cse([Eq(*i) for i in f.items()])
>>> for i in r:
...  print('%s = %s' % i)
...
x0 = sin(t2)
x1 = cos(t2)
x2 = t2_d1**2
x3 = 1/(-t2_d2*x1 + x0*x2)
x4 = t2_d2*x0*x3
x5 = sin(t1)
x6 = cos(t1)
x7 = t1_d1**2
x8 = 1/(-t1_d2*x6 + x5*x7)
x9 = t1_d2*x5*x8
x10 = x1*x2*x3
x11 = x6*x7*x8
x12 = g/(x10 - x11 + x4 - x9)
x13 = x11*x12 + x12*x9
x14 = -x12 - x2_d2

>>> e
[Eq(y1_d2, x13), Eq(y2_d2, x10*x14 + x13 + x14*x4)]

我认为这些例程仅适用于线性变量，因此如果您想求解t1_d1，它将失败。但由于 t1_d1 仅显示为正方形，因此您可以将其替换为 y 并专注于 y。因此，这里有一种关注 t1_d1 和 t2_d2 的方法:

>>> eqs = Tuple(*[Eq(i,0) for i in [eq1, eq2, eq3, eq4, eq16, eq13])
>>> S(focus(eqs.subs(t1_d1**2, y), y, t2_d2)).subs(y, t1_d1**2)
{t1_d1**2: (-L1*t1_d2*sin(t1) + y1_d2)/(L1*cos(t1)),
t2_d2: (-L2*t2_d1**2*cos(t2) - y1_d2 + y2_d2)/(L2*sin(t2))}