Determining the pseudo-time period of a system of $n$-pendulums via Kane's method in Python

We can use Kane's method to integrate the equations of motion for a system of $$n$$ pendulums with arbitrary masses and lengths (see derivation). In particular, if $$(x_i,y_i)$$ denotes the Cartesian coordinates of the $$i$$-th mass $$m_i$$ of length $$l_i$$ of the pendulum for $$1\le i\le n$$, then the Lagrangian formulation is as follows: $$L=T-V=\frac{1}{2}\sum_{i=1}^nm_i(\dot x_i^2+\dot y_i^2)-g\sum_{i=1}^nm_iy_i$$ for $$x_i=\sum_{j=1}^il_j\sin\theta_j$$ and $$y_i=-\sum_{j=1}^il_j\cos\theta_j.$$ When $$n=2$$, this is simply $$L=\frac{1}{2}(m_1+m_2)l_1^2\dot\theta_1^2+\frac{m_2}{2}l_2^2\dot\theta_2^2+m_2l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2)+(m_1+m_2)gl_1\cos\theta_1+m_2gl_2\cos\theta_2,$$ which, up to small angle approximation, is given by $$L=\frac{1}{2}(m_1+m_2)l_1^2\dot\theta_1^2+\frac{m_2}{2}l_2^2\dot\theta_2^2+m_2l_1l_2\dot\theta_1\dot\theta_2-\frac{1}{2}(m_1+m_2)gl_1\theta_1^2+\frac{m_2}{2}gl_2\theta_2^2.$$ Furthermore, the Euler-Lagrange equations $$\frac{\partial L}{\partial\theta_i}-\frac{d}{dt}\left(\frac{\partial L}{\partial \dot\theta_i}\right)=0$$ for $$1\le i\le n$$ yields a system of second-order non-linear differential equations (see here): $$\sum_{j=1}^n\textbf{A}_{ij}\ddot\theta_j=\textbf{b}_i$$ with $$\textbf{A}_{ij}=l_j\cos(\theta_i-\theta_j)\sum_{k=\pi_{ij}}^Nm_k,\\ \textbf{b}_i=-g\sin\theta_i\sum_{j=i}^nm_j-\sum_{j=1}^nl_j\left[\dot\theta_j^2\sin(\theta_i-\theta_j)\right]\sum_{k=\pi_{ij}}^nm_k,\\ \pi_{ij}=\begin{cases}i, &\text{if }i\le j\\ j, & \text{else}\end{cases}.$$ Up to small angle approximation, this is equivalent to $$\sum_{j=1}^nl_j\sum_{k=\pi_{ij}}^Nm_k\ddot\theta_j=-g\theta_i\sum_{j=i}^nm_j-\sum_{j=1}^nl_j\dot\theta_j^2(\theta_i-\theta_j)\sum_{l=\pi_{ij}}^nm_k.$$ In matrix notation, we may write the $$(i,j)$$ entry as: $$\left[l_j\sum_{k=\pi_{ij}}^nm_k\right]\ddot\theta_j=\left[l_j\sum_{k=\pi_{ij}}^nm_k\right]\dot\theta_j^2+\left[gm_j\right]\theta_j.$$ Let $$\bf{\Theta}:=\begin{bmatrix}\theta_1\\ \vdots \\ \theta_n\end{bmatrix},$$ $$\bf M=\begin{bmatrix}l_1\sum_{k=\pi_{11}}^nm_k & & & \bf 0\\ & l_2\sum_{k=\pi_{22}}^nm_k & & & \\ & & \ddots & \\ \bf 0 & & & l_n\sum_{k=\pi_{nn}}^nm_k\end{bmatrix},$$ $$\bf K=\begin{bmatrix}l_1\sum_{k=\pi_{11}}^nm_k & & \bf 0\\ & \ddots & \\ \bf 0 & & l_n\sum_{k=\pi_{nn}}^nm_k\end{bmatrix},$$ and $$\bf L=\begin{bmatrix}gm_1 & & \bf 0\\ & \ddots & \\ \bf 0 & & gm_n\end{bmatrix}.$$ Thus, the equations of motion are $$\textbf{M}\ddot{ \mathbf{\Theta}}=\mathbf{K\dot\Theta}^T\mathbf{\dot\Theta}+\mathbf{L\Theta}.$$ In the case of the $$n$$-pendulum, the solution will consist of oscillations of $$n$$ characteristic frequencies, i.e. normal modes. The normal modes are the real part of the complex-valued vector function $$\mathbf{\Theta}(t)=\begin{bmatrix} \theta_1(t)\\\vdots\\\theta_n(t)\end{bmatrix}=\text{Re}\left(\begin{bmatrix}\mathbf{H}_1\\\vdots\\ \mathbf{H}_n\end{bmatrix}e^{i\omega t}\right)$$ where $$\mathbf{H}_1,\dots\mathbf{H}_n$$ are eigenvectors, $$\omega$$ is the real frequency of the entire system. Up to approximation (I am not sure if this is justified), $${\dot{\mathbf{\Theta}}}^T{\dot{\mathbf{\Theta}}}\approx \mathbf{0}$$ so that the equations of motion become $$\textbf{M}\ddot{\mathbf{\Theta}}=\mathbf{L\Theta}$$, changing the sign of $$\mathbf{L}$$. I made this approximation because I am not certain this system can be linearized using a Jacobian or Taylor development. The value of the normal frequencies are thereby determined as solutons to $$\det(\textbf{K}-\omega^2\textbf{M})=0.$$ That is, $$\det(\textbf{K}-\omega^2\textbf{M})=\prod_{j=1}^n\left[gm_j-\omega^2l_j\sum_{k=\pi_{ij}}^nm_k\right]=0$$ so $$\omega_j\approx \sqrt{\frac{gm_j}{l_j M}}$$ where $$M=\sum_{i=1}^nm_i$$. Therefore, $$\overline{\omega}(n)=\frac{1}{n}\sum_{j=1}^n\sqrt{\frac{gm_j}{l_jM}}$$ so the pseudo-time period as a function of the number of pendulums $$n$$ is approximately $$T(n)=\frac{2\pi n}{\sqrt{\frac{gm_1}{l_1M}}+\dots+\sqrt{\frac{gm_n}{l_nM}}}.$$

Objective: To visualize the relationship between the number of pendulums $$n$$ and the pseudo-time period $$T(n)=\frac{2\pi}{\overline\omega}$$ of the system.

The program implementation below is inspired by jakevdp. However, when I use this approximation method to find the pseudo-time period $$T(n)=\frac{2\pi}{\overline{\omega}}$$, for $$\overline{\omega}=\frac{1}{n}\sum_{i=1}^n\omega_i$$ where $$\omega_i$$ is the characteristic frequency associated with the $$i$$-th mass $$m_i$$ and angle $$\theta_i$$, the graph of T as a function of n seems to be incredibly off in terms of amplitude.

The following implementation defines and solves the equations of motion for a system of n pendulums, with arbitrary masses and lengths (in this case, we let m[i]=1 for all i). This is an initial value problem with theta_1(0)=135,...,theta_n(0)=135 and theta_1'(0)=0,...,theta_n'(0)=0. That is, $$\mathbf{\Theta}(0)=(135,\dots,135)^T$$ and $$\mathbf{\Theta}'(0)=(0,\dots,0)^T.$$

# %matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from sympy import symbols
from sympy.physics import mechanics

from sympy import Dummy, lambdify
from scipy.integrate import odeint

def integrate_pendulum(n, times,
initial_positions=135,
initial_velocities=0,
lengths=None, masses=1):
"""Integrate a multi-pendulum with n sections"""
#-------------------------------------------------
# Step 1: construct the pendulum model

# Generalized coordinates and velocities
# (in this case, angular positions & velocities of each mass)
q = mechanics.dynamicsymbols('q:{0}'.format(n))
u = mechanics.dynamicsymbols('u:{0}'.format(n))

# mass and length
m = symbols('m:{0}'.format(n))
l = symbols('l:{0}'.format(n))

# gravity and time symbols
g, t = symbols('g,t')

#--------------------------------------------------
# Step 2: build the model using Kane's Method

# Create pivot point reference frame
A = mechanics.ReferenceFrame('A')
P = mechanics.Point('P')
P.set_vel(A, 0)

# lists to hold particles, forces, and kinetic ODEs
# for each pendulum in the chain
particles = []
forces = []
kinetic_odes = []

for i in range(n):
# Create a reference frame following the i^th mass
Ai = A.orientnew('A' + str(i), 'Axis', [q[i], A.z])
Ai.set_ang_vel(A, u[i] * A.z)

# Create a point in this reference frame
Pi = P.locatenew('P' + str(i), l[i] * Ai.x)
Pi.v2pt_theory(P, A, Ai)

# Create a new particle of mass m[i] at this point
Pai = mechanics.Particle('Pa' + str(i), Pi, m[i])
particles.append(Pai)

# Set forces & compute kinematic ODE
forces.append((Pi, m[i] * g * A.x))
kinetic_odes.append(q[i].diff(t) - u[i])

P = Pi

# Generate equations of motion
KM = mechanics.KanesMethod(A, q_ind=q, u_ind=u,
kd_eqs=kinetic_odes)
fr, fr_star = KM.kanes_equations(particles, forces)

#-----------------------------------------------------
# Step 3: numerically evaluate equations and integrate

# initial positions and velocities – assumed to be given in degrees

# lengths and masses
if lengths is None:
lengths = np.ones(n) / n
lengths = np.broadcast_to(lengths, n)
masses = np.broadcast_to(masses, n)

# Fixed parameters: gravitational constant, lengths, and masses
parameters = [g] + list(l) + list(m)
parameter_vals = [9.81] + list(lengths) + list(masses)

# define symbols for unknown parameters
unknowns = [Dummy() for i in q + u]
unknown_dict = dict(zip(q + u, unknowns))
kds = KM.kindiffdict()

# substitute unknown symbols for qdot terms
mm_sym = KM.mass_matrix_full.subs(kds).subs(unknown_dict)
fo_sym = KM.forcing_full.subs(kds).subs(unknown_dict)

# create functions for numerical calculation
mm_func = lambdify(unknowns + parameters, mm_sym)
fo_func = lambdify(unknowns + parameters, fo_sym)

# function which computes the derivatives of parameters
def gradient(y, t, args):
vals = np.concatenate((y, args))
sol = np.linalg.solve(mm_func(*vals), fo_func(*vals))
return np.array(sol).T

# ODE integration
return odeint(gradient, y0, times, args=(parameter_vals,))

The above returns generalized physical coordinates, i.e. angular position theta and velocity omega of each pendulum segment, relative to vertical. The following extracts the (x,y) coordinates from the generalized coordinates.

def get_xy_coords(p, lengths=None):
"""Get (x, y) coordinates from generalized coordinates p"""
p = np.atleast_2d(p)
n = p.shape // 2
if lengths is None:
lengths = np.ones(n) / n
zeros = np.zeros(p.shape)[:, None]
x = np.hstack([zeros, lengths * np.sin(p[:, :n])])
y = np.hstack([zeros, -lengths * np.cos(p[:, :n])])
return np.cumsum(x, 1), np.cumsum(y, 1)

Then we fix the number of pendulums to, say, n=10 and determine the pseudo-time period with the following:

n = 10
nperiod = []

#Array containing pseudo-period for a system of n pendulums

for i in range(1,n + 1):
t = np.linspace(0, 10, 1000)
p = integrate_pendulum(i, times=t)

#x, y has first column of zeros

x, y = get_xy_coords(p)
r,s = np.shape(y)

#Call method to find pseudo-period

nperiod.append(computeperiod())

#Takes j, denoting j-th pendulum, as input and returns theta_j for all times where $1\le j\len$. This information corresponds to the j-th column of the y matrix, transformed into polar coordinates.

def theta(j):
theta_j = []
for i in range(0, r):
theta_j.append(math.acos(abs(y[i][j-1]-y[i][j])))
#We should technically divide by the length of the pendulum in abs(.)

timenew = [i for i in range(1,r + 1)]
graph_j = pd.Series(data=theta_j, index=timenew)

#Returns array omega_j which is the numerical time derivative of all values contained in the r x 1 array theta_j, i.e. the value of the j-th angle theta_j for all times

#Returns pseudo-time period using formula T=2π/omega

def computeperiod():
series = []
for j in range(1,s):
series.append(2 * math.pi/(theta(j).mean()))
numberline = [i for i in range(1,s)] #Here s=n+1
timeperiod = pd.Series(data=series, index=numberline)
return abs(timeperiod.mean())
print(nperiod)
numbers = [i for i in range(len(nperiod))]
finalperiods = pd.Series(data=nperiod, index=numbers)
finalperiods.plot()

When we plot it, we get the following graph for n=10: However, it should behave asymptotically according to $$T(n)\approx 2\pi n^{3/2}\sqrt{l/g}$$, where $$l = 1$$, as below. Furthermore, its amplitude should be reduced by a factor of approximately $$100,000$$. Thus, it seems as though something is clearly wrong with this formulation. I suppose there is the possibility that all of the small angle approximations produced this discrepancy; however, it seems to be an error in physical calculation. Note, the run time for this implementation is also exponential, which doesn't help with visualization.

Is the pseudo-period of an $$n$$-body system of pendulums given by $$T(n)=\frac{2\pi n}{\sqrt{\frac{gm_1}{l_1M}}+\dots+\sqrt{\frac{gm_n}{l_nM}}}?$$ That is, is the above derivation correct and, if so, what is the reason for the large discrepancy between the expected result and the simulation?

Any help would be much appreciated!