Consider the pendulum suspended on the pendulum. We weill use dAlembert principle to derive eqaution of motion in generalized coordinaes. Naturally, we choose two angles as coordinates which comply with contraints.
var('t')
var('l1 l2 m1 m2 g')
xy_names = [('x1','x1'),('y1','y1'),('x2','x2'),('y2','y2')]
uv_names = [('phi1','\\varphi_1'),('phi2','\\varphi_2')]
load('cas_utils.sage')
to_fun, to_var = make_symbols(xy_names,uv_names)
x2u = {x1:l1*sin(phi1),\
y1:-l1*cos(phi1),\
x2:l1*sin(phi1)+l2*sin(phi2),\
y2:-l1*cos(phi1)-l2*cos(phi2)}
transform_virtual_displacements(xy_names,uv_names,verbose=True)
dAlemb = (m1*x1.subs(x2u).subs(to_fun).diff(t,2))*dx1_polar + \
(m1*y1.subs(x2u).subs(to_fun).diff(t,2)+m1*g)*dy1_polar+\
(m2*x2.subs(x2u).subs(to_fun).diff(t,2))*dx2_polar + \
(m2*y2.subs(x2u).subs(to_fun).diff(t,2)+m2*g)*dy2_polar
dAlemb = dAlemb.subs(to_var)
showmath(dAlemb)
eq1 = dAlemb.expand().coefficient(dphi1).trig_simplify()
eq2 = dAlemb.expand().coefficient(dphi2).trig_simplify()
showmath(eq1)
showmath(eq2)
sol = solve([eq1,eq2],[phi1dd,phi2dd])[0]
showmath(sol)
showmath( sol[0].rhs().denominator() )
(l1*(2*m1+m2-m2*cos(2*phi1-2*phi2))).expand_trig().expand_trig().expand().show()
bool ( -2*sol[0].rhs().denominator()==(l1*(2*m1+m2-m2*cos(2*phi1-2*phi2))).expand_trig().expand_trig().expand() )
Since the "textbook" solution contains a slightly different form, let's check if we have these formulas:
$$T(\varphi_1,\varphi_2,\dot{\varphi}_1,\dot{\varphi}_2) = \frac{m_1}{2} l_1^2 \dot{\varphi}_1^2 + \frac{m_2}{2} \left( l_1^2 \dot{\varphi}_1^2 + l_2^2 \dot{\varphi}_2^2 + 2 l_1 l_2 \dot{\varphi}_1 \dot{\varphi}_2 \cos(\varphi_1-\varphi_2) \right)$$$$V(\varphi_1,\varphi_2) = -(m_1+m_2) g l_1 \cos\varphi_1 - m_2 g l_2 \cos\varphi_2$$ $$m_{2}l_{2}\ddot{\varphi}_{2}\cos\left(\varphi_{1}-\varphi_{2}\right)+\left(m_{1}+m_{2}\right)l_{1}\ddot{\varphi}_{1}+m_{2}l_{2}\dot{\varphi}_{2}^{2}\sin\left(\varphi_{1}-\varphi_{2}\right)+\left(m_{1}+m_{2}\right)g\sin\varphi_{1}=0$$ $$l_{2}\ddot{\varphi}_{2}+l_{1}\ddot{\varphi}_{1}\cos\left(\varphi_{1}-\varphi_{2}\right)-l_{1}\dot{\varphi}_{1}^{2}\sin\left(\varphi_{1}-\varphi_{2}\right)+g\sin\varphi_{2}=0$$
rown_wiki = [m2*l2*cos(phi1-phi2)*phi2dd+(m1+m2)*l1*phi1dd+m2*l2*phi2d^2 * sin(phi1-phi2)+ (m1+m2)*g*sin(phi1),\
l2*phi2dd+l1*cos(phi1-phi2)*phi1dd-l1*phi1d^2*sin(phi1-phi2)+g*sin(phi2)]
showmath(rown_wiki[0])
showmath(rown_wiki[1])
rown_wiki[0].show()
(eq1/l1).reduce_trig().show()
rown_wiki[0].show()
bool((eq1/l1) == rown_wiki[0] )
(eq2/l2/m2).reduce_trig().show()
rown_wiki[1].show()
bool((eq2/l2/m2) == rown_wiki[1] )
Ekin = 1/2*(m1*x1.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2+\
m1*y1.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2+\
m2*x2.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2+\
m2*y2.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 )
Epot = m1*g*y1.subs(x2u)+m2*g*y2.subs(x2u)
showmath( Epot.collect(cos(phi1)) )
showmath( Epot )
showmath( Ekin.trig_simplify() )
L = Ekin - Epot
len(L.expand().operands())
EL1 = L.diff(phi1d).subs(to_fun).diff(t).subs(to_var) - L.diff(phi1)
EL2 = L.diff(phi2d).subs(to_fun).diff(t).subs(to_var) - L.diff(phi2)
EL1.expand().operands()
EL1 = (EL1/l1).trig_reduce()
EL2 = (EL2/l2).trig_reduce()
showmath(EL1)
sol = solve([EL1,EL2],[phi1dd,phi2dd])[0]
show(sol)
expr = sol[0].rhs()
for ex_ in expr.factor().numerator().operands():
show(ex_)
import numpy as np
ode = [phi1d,phi2d]+[sol[0].rhs(),sol[1].rhs()]
ode = map(lambda x:x.subs({l1:1,l2:1,m1:1,m2:1,g:9.81}),ode)
times = srange(0,5,.01)
numsol = desolve_odeint(ode,[2.1,0,0,0],times,[phi1,phi2,phi1d,phi2d])
p = line ( zip(np.sin(numsol[:,0])+np.sin(numsol[:,1]),\
-np.cos(numsol[:,0])-np.cos(numsol[:,1])), color='gray' )
p.show(figsize=4)
def plot_dp(f1,f2,pars):
mass1 = vector([x1,y1]).subs(x2u).subs(pars).subs({phi1:f1,phi2:f2})
mass2 = vector([x2,y2]).subs(x2u).subs(pars).subs({phi1:f1,phi2:f2})
plt = point([(0,0),mass1],aspect_ratio=1,size=40)
plt += point(mass2,xmin=-2,xmax=2,ymin=-2,ymax=2,size=40)
plt += line([(0,0),mass1,mass2],color='gray')
return plt
plot_dp(numsol[213,0],numsol[213,1],{l1:1,l2:1})+p
@interact
def _(ith=slider(0,numsol.shape[0]-1)):
f1,f2 = numsol[ith,:2]
plot_dp(f1,f2,{l1:1,l2:1}).show(axes=False)
\newpage