sympy：双摆方程组

Question

我想求解文章 Double Pendulum 中描述的方程组，与作者相似，在最后一步他声称他 "use a computer algebra program to solve equations (13) and (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())

正在研究类似如何在 SymPy 中求解线性方程组的问题？，我预计 SymPy 会简化每个符号的 return 方程，但是对于方程 (16)，我得到的结果仅为：L1、g、t1、L2、t2，而不是对于 ~~y2_d2~~ t2_d2。对于等式 (13)，我得到了例外。

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

已编辑： 对于 solve/find 两个未知数 ~~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))}

Answer 2

多亏了 smichr 的线索，我才能够得到想要的结果。

from sympy import symbols, sin, cos, solve, simplify, Eq

m1, m2 = symbols('m1, m2')
g = symbols('g')
t1, t2, t1_d1, t1_d2, t2_d1, t2_d2 = symbols('t1, t2, t1_d1, t1_d2, t2_d1, t2_d2')
L1, L2 = symbols('L1, L2')

x1_d2 = -t1_d1**2*L1*sin(t1) + t1_d2*L1*cos(t1)
y1_d2 =  t1_d1**2*L1*cos(t1) + t1_d2*L1*sin(t1)
x2_d2 =  x1_d2 - t2_d1**2*L2*sin(t2) + t2_d2*L2*cos(t2)
y2_d2 =  y1_d2 + t2_d1**2*L2*cos(t2) + t2_d2*L2*sin(t2)
eq13 = Eq(sin(t1)*(m1*y1_d2 + m2*y2_d2 + m2*g + m1*g) + cos(t1)*(m1*x1_d2 + m2*x2_d2), 0)
eq16 = Eq(sin(t2)*(m2*y2_d2 + m2*g) + cos(t2)*(m2*x2_d2), 0)


solve([eq13, eq16], [t1_d1, t2_d2], dict=True)

结果：

[{t1_d1: -sqrt((L1*m2*t1_d2*sin(2*t1 - 2*t2) - (2*L1*m1*t1_d2 + L1*m2*t1_d2 + 2*L2*m2*t2_d1**2*sin(t1 - t2) + 2*g*m1*sin(t1) + g*m2*sin(t1) + g*m2*sin(t1 - 2*t2))*tan(2*t1 - 2*t2))/(L1*m2*sin(2*t1 - 2*t2)*tan(2*t1 - 2*t2))),
  t2_d2: -(L1*m1*t1_d2 + L1*m2*t1_d2 + L2*m2*t2_d1**2*sin(t1 - t2) + g*m1*sin(t1) + g*m2*sin(t1))/(L2*m2*cos(t1 - t2))},
 {t1_d1: sqrt((L1*m2*t1_d2*sin(2*t1 - 2*t2) - (2*L1*m1*t1_d2 + L1*m2*t1_d2 + 2*L2*m2*t2_d1**2*sin(t1 - t2) + 2*g*m1*sin(t1) + g*m2*sin(t1) + g*m2*sin(t1 - 2*t2))*tan(2*t1 - 2*t2))/(L1*m2*sin(2*t1 - 2*t2)*tan(2*t1 - 2*t2))),
  t2_d2: -(L1*m1*t1_d2 + L1*m2*t1_d2 + L2*m2*t2_d1**2*sin(t1 - t2) + g*m1*sin(t1) + g*m2*sin(t1))/(L2*m2*cos(t1 - t2))}]

sympy：双摆方程组

sympy: system of equations for double pendulum

python

linear-algebra

solver

sympy