$N$ pendula¶

import numpy as np 

var('t,g')

load('cas_utils.sage')

N = 3
var('t g')
for i in range(1,1+N):
    var('l%d m%d'%(i,i))

xy_names = [('x%d'%i,'x_%d'%i) for i in range(1,1+N)] 
xy_names += [('y%d'%i,'y_%d'%i) for i in range(1,1+N)] 
uv_names = [ ('phi%d'%i,'\\varphi_%d'%i) for i in range(1,1+N)] 

load('cas_utils.sage')

to_fun, to_var = make_symbols(xy_names,uv_names)

ls = [vars()['l%d'%i] for i in range(1,1+N)]
xs = [vars()['x%d'%i] for i in range(1,1+N)]
ys = [vars()['y%d'%i] for i in range(1,1+N)]
ms = [vars()['m%d'%i] for i in range(1,1+N)]
phis = [vars()['phi%d'%i] for i in range(1,1+N)]
phids = [vars()['phi%dd'%i] for i in range(1,1+N) ]
phidds = [vars()['phi%ddd'%i] for i in range(1,1+N) ]

showmath(phis)

x2u = {x1:l1*sin(phi1),\
       y1:-l1*cos(phi1) }

for x_prev,x_,y_prev,y_,l_,phi_ in zip(xs[:-1],xs[1:],ys[:-1],ys[1:],ls[1:],phis[1:]):
    x2u[x_] = x2u[x_prev] + l_*sin(phi_)
    x2u[y_] = x2u[y_prev] - l_*cos(phi_)

transform_virtual_displacements(xy_names,uv_names,verbose=False)

dxs = [vars()['dx%d_polar'%i] for i in range(1,1+N) ]
dys = [vars()['dy%d_polar'%i] for i in range(1,1+N) ]

dAlemb = sum( (m_*x_.subs(x2u).subs(to_fun).diff(t,2))*dx_ for m_,x_,dx_ in zip(ms,xs,dxs) )
dAlemb += sum( (m_*x_.subs(x2u).subs(to_fun).diff(t,2) + m_*g)*dx_ for m_,x_,dx_ in zip(ms,ys,dys) )
dAlemb = dAlemb.subs(to_var)

#showmath(dAlemb)

dphis = [vars()['dphi%d'%i] for i in range(1,1+N) ]

eqs = [dAlemb.expand().coefficient(dphi_).trig_simplify() for dphi_ in dphis]

showmath(eqs[1].trig_reduce())

sol =  solve(eqs,phidds)[0]

len(sol[1].rhs().trig_reduce().operands())

pars= {m_:1 for m_ in ms}

for i,l_ in enumerate(ls):
    pars[l_]  = 1/ (i+1)
pars[g] = 9.81

phidds

ode =  phids + [sol_.rhs() for sol_ in sol]
ode = map(lambda x:x.subs(pars),ode)

times = srange(0,5,.01)

ics = [0]*(N*2)
ics[-1] = 33.01
ics[:N]= [0*pi.n()]*N
numsol = desolve_odeint(ode,ics,times, phis + phids)

phi_subs = lambda ith: {phi_:numval_ for phi_,numval_ in zip(phis,numsol[ith,:N])}

@interact
def _(ith=slider(0,numsol.shape[0]-1,step_size=10)):
    xnum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in xs]
    ynum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in ys]
    plt = line(zip(xnum,ynum),\
        xmin=-N,xmax=N,ymin=-N,ymax=N,marker='o')
    plt.show(axes=False,figsize=5,aspect_ratio=1)

Ekin = 1/2 * sum( 
        m_*x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 +\
        m_*y_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 \
                 for m_,x_,y_ in zip(ms,xs,ys))

Ekin = Ekin.trig_simplify()
Epot = sum(m_*g*y_.subs(x2u) for m_,y_ in zip(ms,ys))

showmath(Ekin)

showmath(Epot)

L = Ekin - Epot

ELs = [L.diff(phid_).subs(to_fun).diff(t).subs(to_var) - L.diff(phi_)\
       for (phi_,phid_) in zip(phis,phids)]

ELs = [(EL_/l_).trig_reduce() for EL_,l_ in zip(ELs,ls)]

showmath(ELs[0])

sol =  solve(ELs,phidds)[0]
#show(sol)

pars= {m_:1 for m_ in ms}

for i,l_ in enumerate(ls):
    pars[l_]  = 0.3+1/ (i+1)
pars[g] = 9.81

ode =  phids + [sol_.rhs() for sol_ in sol]
ode = map(lambda x:x.subs(pars),ode)

times = srange(0,5,.01)

ics = [0]*(N*2)
ics[-1] = 33.01
ics[:N]= [0*pi.n()]*N
numsol = desolve_odeint(ode,ics,times, phis + phids)

phi_subs = lambda ith: {phi_:numval_ for phi_,numval_ in zip(phis,numsol[ith,:N])}

@interact
def _(ith=slider(0,numsol.shape[0]-1,step_size=10)):
    xnum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in xs]
    ynum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in ys]
    plt = line(zip(xnum,ynum),\
        xmin=-N,xmax=N,ymin=-N,ymax=N,marker='o')
    plt.show(axes=False,figsize=5,aspect_ratio=1)