sympy: system of equations for double pendulum

Question

I want to solve system of equation described in the article Double Pendulum, similar like author, where in last step he claims that he "use a computer algebra program to solve equations (13) and (16)". To do this, I rewrite this equations in 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())

Looking on similar question like How can I solve system of linear equations in SymPy? , I expected that SymPy simplify, and return equations in respect to each symbol, however for equation (16) I got the result only for: L1, g, t1, L2, t2, and not for ~~y2_d2~~ t2_d2. And for equation (13) I got 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

Edited: How should my code look like, to solve/find two unknowns ~~y1_d2, y2_d2~~ t1_d1, t2_d2 - equations (13, 16)?

smichr · Answer 1 · 2019-08-13T14:06:29.337

This sort of manipulation is described here; using the focus routine described there you can feed your full equation set (cast as Eq instances) to focus and specify y1_d2 and y2_d2 as the variables on which you want to focus:

>>> 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

A compact result can be obtained with 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)]

I think the routines will only work with linear variables so if you want to solve for t1_d1 it will fail. But since t1_d1 only appears as a square you can replace it with y and focus on y. So here is a way to focus on t1_d1 and 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))}

I made mistake in my question. Instead finding y1_d2, y2_d2 I need t1_d1, t2_d2, so I rearrange your code `F = (t1_d2, t2_d2); focus([Eq(i, 0) for i in [eq1, eq2, eq3, eq4, eq16, eq13]], *F)` and I got exception: `'BooleanFalse' object has no attribute 'lhs'`, but when I swap `eq13` and `eq16` no exception occur. Does order of formulas matters? — NiegodziwyBeru, Aug 13 '19 at 08:44
It shouldn't matter. I updated the routines to watch out for this situation. Again, lightly tested (and it's been a while since I wrote those). Does the updated result look like what you were looking for? — smichr, Aug 13 '19 at 14:07
Not quit, because e.g. for `t2_d2` I hope that SymPy get rid of `y2_d2`, `y1_d1`. Basically I was expected that every previous equation will be rearranged, put to next one and at the end solved. So I end up with constants (g, m, L) and theta with his derivatives. But your comment gave my a clue - *subs* and *Eq*, and after some reading I was able to rewrite my code. Thanks. — NiegodziwyBeru, Aug 15 '19 at 17:26

score 2 · Accepted Answer · answered Aug 15 '19 at 17:29

Thanks to smichr clue I was able to get desired result.

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)

Result:

[{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: system of equations for double pendulum

2 Answers2