0
$\begingroup$

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
y0 = np.deg2rad(np.concatenate([np.broadcast_to(initial_positions, n),
                                np.broadcast_to(initial_velocities, n)]))

# 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[0]

# 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[1] // 2
    if lengths is None:
        lengths = np.ones(n) / n
    zeros = np.zeros(p.shape[0])[:, 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

    return pd.Series(data=np.gradient(graph_j.values), index=graph_j.index) 

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

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.