Automatically Deriving the Angular Accelerations for N-Pendula

n-pendula.py

from sympy import *
import matplotlib.pyplot as plt

# time variable
t = symbols('t')

# number of pendula + 1
Z = 3 + 1

# Constants
G = symbols('G')
M = symbols('M0:%s' % Z)
L = symbols('L0:%s' % Z)

# Angles and their first and second derivatives
u = symbols('u0:%s' % Z)
v = symbols('v0:%s' % Z)
w = symbols('w0:%s' % Z)

# cartesian coordinates
x, y = [0], [0]
for i in xrange(1,Z):
    x.append( x[-1] + L[i] * cos( u[i](t) ) )
    y.append( y[-1] + L[i] * sin( u[i](t) ) )

# velocity in cartesian coordinates
K = 0
for i in xrange(1,Z):
    K += M[i] * (diff(x[i], t) ** 2 + diff(y[i], t) ** 2 ) / 2
    
K = simplify( K )

P = 0
for i in xrange(1,Z):
    P += M[i] * G * y[i]

L = K - P
print L

EL = [0]
for i in xrange(1,Z):
    EL.append( L.diff( Derivative(u[i](t), t) ).diff(t) - L.diff(u[i](t) ) )

S = {}
for i in xrange(1,Z):
    S[Function(u[i].name)(t)] = u[i]
    S[Derivative(Function(u[i].name)(t), t)] = v[i]
    S[Derivative(Function(u[i].name)(t), t, t)] = w[i]

for i in xrange(1,Z):    
    EL[i] = EL[i].subs(S)

first_order_odes = solve(EL[1:], w[1:])
Q = [0]
for i in xrange(1,Z):
    Q.append( simplify(first_order_odes[w[i]]) )
    
things_to_replace = 'uvML'
for thing in things_to_replace:
    for i in xrange(1,Z):
        for j in xrange(1,Z):
            Q[j] = str(Q[j]).replace('%s%s' % (thing, i), '%s[%s]' % (thing, i))
        
for i in xrange(1,Z):
    print 'Q%s' % i
    print Q[i]

Angular Accelerations for Triple Pendulum (i.e. d_tt theta_i(t)). u == theta_i, v = d/dt u, w = d/dt v .

Angular accelerations
Q1
(M[3]*(M[2] + M[3])*(cos(u[1] - u[3]) - cos(u[1] - 2*u[2] + u[3]))*(-G*cos(u[3]) + L[1]*v[1]**2*sin(u[1] - u[3]) + L[2]*v[2]**2*sin(u[2] - u[3]))/2 + (M[3]*cos(u[1] - u[3])*cos(u[2] - u[3]) + (-M[2] - M[3])*cos(u[1] - u[2]))*(G*M[2]*cos(u[2]) + G*M[3]*cos(u[2]) - L[1]*M[2]*v[1]**2*sin(u[1] - u[2]) - L[1]*M[3]*v[1]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[2] - u[3])) + (M[2] - M[3]*cos(u[2] - u[3])**2 + M[3])*(G*M[1]*cos(u[1]) + G*M[2]*cos(u[1]) + G*M[3]*cos(u[1]) + L[2]*M[2]*v[2]**2*sin(u[1] - u[2]) + L[2]*M[3]*v[2]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[1] - u[3])))/(L[1]*(-2*M[3]*(M[2] + M[3])*cos(u[1] - u[2])*cos(u[1] - u[3])*cos(u[2] - u[3]) + M[3]*(M[2] + M[3])*cos(u[1] - u[3])**2 + M[3]*(M[1] + M[2] + M[3])*cos(u[2] - u[3])**2 + (M[2] + M[3])**2*cos(u[1] - u[2])**2 - (M[2] + M[3])*(M[1] + M[2] + M[3])))
Q2
(-M[3]*((M[2] + M[3])*cos(u[1] - u[2])*cos(u[1] - u[3]) - (M[1] + M[2] + M[3])*cos(u[2] - u[3]))*(-G*cos(u[3]) + L[1]*v[1]**2*sin(u[1] - u[3]) + L[2]*v[2]**2*sin(u[2] - u[3])) + (M[3]*cos(u[1] - u[3])*cos(u[2] - u[3]) - (M[2] + M[3])*cos(u[1] - u[2]))*(G*M[1]*cos(u[1]) + G*M[2]*cos(u[1]) + G*M[3]*cos(u[1]) + L[2]*M[2]*v[2]**2*sin(u[1] - u[2]) + L[2]*M[3]*v[2]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[1] - u[3])) + (M[1] + M[2] - M[3]*cos(u[1] - u[3])**2 + M[3])*(G*M[2]*cos(u[2]) + G*M[3]*cos(u[2]) - L[1]*M[2]*v[1]**2*sin(u[1] - u[2]) - L[1]*M[3]*v[1]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[2] - u[3])))/(L[2]*(-2*M[3]*(M[2] + M[3])*cos(u[1] - u[2])*cos(u[1] - u[3])*cos(u[2] - u[3]) + M[3]*(M[2] + M[3])*cos(u[1] - u[3])**2 + M[3]*(M[1] + M[2] + M[3])*cos(u[2] - u[3])**2 + (M[2] + M[3])**2*cos(u[1] - u[2])**2 - (M[2] + M[3])*(M[1] + M[2] + M[3])))
Q3
(-(M[2] + M[3])*(cos(u[1] - u[3]) - cos(u[1] - 2*u[2] + u[3]))*(G*M[1]*cos(u[1]) + G*M[2]*cos(u[1]) + G*M[3]*cos(u[1]) + L[2]*M[2]*v[2]**2*sin(u[1] - u[2]) + L[2]*M[3]*v[2]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[1] - u[3])) - 2*(M[2] + M[3])*(-G*cos(u[3]) + L[1]*v[1]**2*sin(u[1] - u[3]) + L[2]*v[2]**2*sin(u[2] - u[3]))*(M[1] + M[2] + M[3] - (M[2] + M[3])*cos(u[1] - u[2])**2) + 2*((M[2] + M[3])*cos(u[1] - u[2])*cos(u[1] - u[3]) - (M[1] + M[2] + M[3])*cos(u[2] - u[3]))*(G*M[2]*cos(u[2]) + G*M[3]*cos(u[2]) - L[1]*M[2]*v[1]**2*sin(u[1] - u[2]) - L[1]*M[3]*v[1]**2*sin(u[1] - u[2]) + L[3]*M[3]*v[3]**2*sin(u[2] - u[3])))/(2*L[3]*(-2*M[3]*(M[2] + M[3])*cos(u[1] - u[2])*cos(u[1] - u[3])*cos(u[2] - u[3]) + M[3]*(M[2] + M[3])*cos(u[1] - u[3])**2 + M[3]*(M[1] + M[2] + M[3])*cos(u[2] - u[3])**2 + (M[2] + M[3])**2*cos(u[1] - u[2])**2 - (M[2] + M[3])*(M[1] + M[2] + M[3])))