Elastic Pendulum in Python

Elastic Pendulum is a simple mechanical structure, which is composed of only one spring and a mass. Unlike its simple structure, it also has interesting and complex motion and meaning. For example, in addition to the cyclic motion trajectory, it also has chaotic motion. On the other hand, this model is very similar to the SLIP (Spring-Loaded Inverted Pendulum) model [1] commonly seen in robotics.

Before we start…

In this article, the construction of the elastic pendulum model will be described in detail. The programs used to simulate the elastic pendulum in this article is available on my GitHub. If you need them, please click the link below:

GitHub - SheenaEita/Elastic_Pendulum

In addition to the Python version, the above link also includes a MATLAB version. I personally think that the results drawn by MATLAB look better, but there is no difference in what they express. Therefore, this article will also use some pictures or animations generated by MATLAB for illustration.

SLIP model and the pendulum are similar. If you are interested, here are the articles I’ve written about the SLIP model in the past.

Construct Spring-Loaded Inverted Pendulum(SLIP) Model in MATLAB

If you have any other related questions or needs, please feel free to contact me :)

Email : aithereita@gmail.com

Content:
1. Package, Parameters and Coordinates
2. Lagrangian and Equation of Motion
3. Simulation and Analysis
4. Conclution and Prospect

1. Package, Parameters and Coordinates

To build a simulation model, symbolic operations will be needed, so import “sympy” [2] at the beginning of the program to help us complete the equation derivation. Similarly, to solve the equation, ODE integration instructions are needed. In order to make the results more intuitive and understandable, we also need some drawing tools to visualize the data.

import numpy as np
import sympy as sym
from scipy.integrate import odeint as inte
import matplotlib.pyplot as plt

import math
from math import pi

import csv

For an elastic pendulum system, there are only few parameters: time t, mass m, gravity constant g, elastic constant of the spring k, and the spring’s initial length l.

t, m, g, k, l = sym.symbols('t, m, g, k, l')
r, theta = sym.symbols(r'r \theta', cls=sym.Function)

I define the coordinates of this system to be the same as the standard polar coordinates [r, θ]. The fixed point of the spring is at the origin of the coordinates [0, 0]. Next, describing r and θ as functions of time, and defining their first and second derivatives. These will be used in the derivation of equation of motion.

r = r(t)
theta = theta(t)
r_d = sym.diff(r, t)
theta_d = sym.diff(theta, t)
r_dd = sym.diff(r_d, t)
theta_dd = sym.diff(theta_d, t)

I let the program output part of the results. The result of θ double dot is:

Similarly, in order to make it easier for us to understand, define the Cartesian coordinate [x, y], and establish the conversion between it and the polar coordinates. Through this conversion, we can convert the system in the polar coordinate system to the Cartesian coordinate system for simulation and drawing.

x = r * sym.cos(theta)
y = r * sym.sin(theta)

x_d = sym.diff(x, t)
y_d = sym.diff(y, t)

2. Lagrangian and Equation of Motion

For such a complex motion system, Lagrangian Mechanism is usually used to calculate its motion trajectory. First, we describe the kinetic and potential energy of the system, and then subtract them to get the Lagrangian.

T = 0.5*m*(x_d**2 + y_d**2);
# T = 1/2 * m * (x_dot^2 + y_dot^2);

V = m*g*y + 0.5*k*(r-l)**2;
# V = mgz + 1/2 * k * (r-l)^2 ;

L = T - V;

Lagrangian Equation

LE1 = sym.diff(L, r) - sym.diff(sym.diff(L, r_d), t).simplify()
LE2 = sym.diff(L, theta) - sym.diff(sym.diff(L, theta_d), t).simplify()

Python also provides the results of symbolic operations. The results are as follows:

This is really a very long equation. I am so glad that I was born in an era where using computers is so convenient.

In order to let the computer to do numerical calculations, the above equations need to be organized into the state space form.

EP_sim = sym.solve([LE1, LE2], (r_dd, theta_dd))

d2rdt2_f = sym.lambdify((t, m, g, k, l, r, theta, r_d, theta_d), EP_sim[r_dd])
d2thetadt2_f = sym.lambdify((t, m, g, k, l, r, theta, r_d, theta_d), EP_sim[theta_dd])
drdt_f = sym.lambdify(r_d, r_d)
dthetadt_f = sym.lambdify(theta_d, theta_d)

def dSdt(S, t, m, g, k, l):
    r, theta, r_d, theta_d = S
    return [
        drdt_f(r_d),
        dthetadt_f(theta_d),
        d2rdt2_f(t, m, g, k, l, r, theta, r_d, theta_d),
        d2thetadt2_f(t, m, g, k, l, r, theta, r_d, theta_d),
    ]

3. Simulation and Analysis

We have completed the model building, and then, as long as the values of parameters and the initial conditions of motion are set, the corresponding motion trajectory can be obtained through numerical calculation of the computer.

# parameters of the system
t = np.linspace(0, 30, 10001) # s
g = 9.81 #m/s^2
k = 4*math.pi**2 # N/m
m = 1 # kg
l = 1 # m

#initial condition
ic_r = 1
ic_theta = -60*pi/180
ic_r_d = 0
ic_theta_d = 0

ans = inte(dSdt, y0=[ic_r, ic_theta, ic_r_d, ic_theta_d], t=t, args=(m, g, k, l))

r — t plot

θ — t plot

The results obtained by program are expressed in polar coordinates. In order to draw the pendulum’s motion, we need to convert them into Cartesian coordinates.

x_result = np.empty((10001))
y_result = np.empty((10001))
for i in range(1, 10001):
  x_result[i] = r_result[i]*sym.cos(theta_result[i])
  y_result[i] = r_result[i]*sym.sin(theta_result[i])

After the above steps, we have completed the whole process from definition to display of an elastic pendulum system from polar coordinates to Cartesian coordinates and from mathematical expression to visualization. Finally, save the calculation results as a csv file to retain the data for future use.

plt.plot(t, r_result)
plt.plot(t, theta_result)
plt.plot(x_result[1:10000], y_result[1:10000])

file_name = 'elastic_pendulum_sim_result.csv'
with open(file_name, 'w', newline='') as csvFile:
  csvWriter = csv.writer(csvFile)
  csvWriter.writerow(t)
  csvWriter.writerow(r_result)
  csvWriter.writerow(theta_result)
  csvWriter.writerow(x_result)
  csvWriter.writerow(y_result)

4. Conclution and Prospect

The pendulum is a simple mechanical system, which is also a good example to practice Lagrangian mechanism and explore chaotic systems. There are lots of knowledge about this on the internet. By changing the elastic coefficient of the spring, setting different weights, or altering the initial position or velocity, you can get completely different trajectory.

If we use programs to construct a simulated model in the computer, it will help us understand the motion of this system under different conditions more conveniently. I have built a physical pendulum mechanism and added an LED at the end in order to obtain the motion trajectory by long exposure. The result was that I managed to see the motion trajectory, but compared to computer simulation, I was unable to set the initial conditions and prepare ideal experimental environment accurately.

Building a 2D single pendulum model was not enough for me, so I also provided the 3D version I built using MATLAB on my GitHub. If you are interested, you can visit my GitHub to get the related programs. Here is the link again:

GitHub - SheenaEita/Elastic_Pendulum

And for python version, I wrote it on Google colab, and the link is as below:

Google Colaboratory

Thank you for reading. I’m honored to have had the opportunity to share with you.

Reference

[1] R. M. Alexander, Elastic mechanisms in animal movement. Cambridge University Press, 1988.

[2] Meurer, Aaron, et al. SymPy: symbolic computing in Python. PeerJ Computer Science 3, 2017.