Simulating a pendulum using Python

Aim

The aim of this program is to solve the ordinary differential equation(ODE) which governs the motion of a pendulum and to simulate its motion.

Equation

The ODE equation is `(d^2theta)/dt^2 + b/m*(d theta)/dt+g/L*sintheta=0`which is split into two first order differential equations `(d theta_1)/dt=theta_2` and `(d theta_2)/dt=-b/m*theta_2-g/L*sin theta_1`

Program

In the program the length of the pendulum(l) is 1m, the mass of the ball(m) is 1kg, damping coefficent(b) is 0.05 and acceleration due to gravity is taken as 9.81 m/s^2. The initial values of angular displacement and angular velocity is taken as 0 and 3 rad/s. The time over which the ODE is integrated is 20s.

'odeint' from scipy was used to integrate the system of differential equations. The function 'model' was defined to take the pendulum parameters and inital conditions and returned the differential equation for 'odeint' to use. Angular displacement and angular velocity values were stored in 'theta' array. The variation of angular displacement and angular velocity with time are plotted. The angular displacement values are used to plot the moving end of the string. A circle is also plotted with its center at a distance of 'radius' from the moving end of the string where 'radius' is the radius of the circle. Each plot is saved using 'savefig' command and is used to make a gif.

Results

Code

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#function whuch returns the differential equation
def model(theta_init,t,b,g,l,m):
	theta1=theta_init[0]
	theta2=theta_init[1]
	dtheta1_dt=theta2
	dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return dtheta_dt
#pendulum parameters
l=1
b=0.05
g=9.81
m=1

theta_init=[0,3]#initial conditions
t=np.linspace(0,20,150)

angles = np.linspace(0, 2*math.pi, 100) #angles to plot the circle
radius=0.1
count=1;


theta=odeint(model,theta_init,t,args=(b,g,l,m))
ang_disp=theta[:,0] #stores the angular displacement values

#plotting the variation of angular displacement and velocity with time
plt.figure()
plt.plot(t,theta[:,0],color="red")
plt.plot(t,theta[:,1],color="blue")
plt.xlabel("Time")
plt.title("Variation of angular displacement and velocity with time")
plt.legend(["Angular displacement","Angular velocity"])
plt.show()

#25 images are created for the gif
for i in range(25):
	plt.figure()
	#plotting the string and the moving end
	x_coord=l*math.sin(ang_disp[i])
	y_coord=3-l*math.cos(ang_disp[i])
	plt.plot([0,x_coord],[3,y_coord])

	#plotting the circle
	x_circle = x_coord+radius*math.sin(ang_disp[i])+radius*np.cos(angles)
	y_circle = y_coord-radius*math.cos(ang_disp[i])+radius*(np.sin(angles))
	plt.plot(x_circle,y_circle)

	plt.xlim(-2,2)
	plt.ylim(0,3)
	plt.savefig(f"pendulum{count}") #saving each figure
	count=count+1