Solving the ODE for simple pendulum followed by its simulation

Objective:

• To solve the following ODE of a simple pendulum with damping and then plot the displacement and velocity values on a graph

`(d^2theta)/dt^2 + b/m(d theta)/dt + g/lsintheta=0`

where

g = gravity in `m/s^2`,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

• And then simulate the simple pendulum up to 20 seconds

Problem Statement:

How does displacement and velocity vary with time and how can these values be used to simulate the pendulum

Solution Procedure:

As we want the angular displacement and angular velocity values stating from the initial position of 0 degrees and 3radians respectively, we have initialised those conditions also. And finally using the linspace command, created an array of 200 values between 0 and 20 corresponding to each time element we will be plotting the values in.

Now, by the use of the integration function from the scipy.integrate module, I have calculated each and every value of angular displacement and angular velocity up to 20 seconds and stored them in 2 columns of an array

Lastly, to complete the 1^st part of the objective, I have used the plotting commands to plot the angular displacement vs time and angular velocity vs time graphs

Here is the detailed code for this:

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


#function that returns dz/dt
def model(theta,t,b,m,g,l):
	theta1=theta[0]
	theta2=theta[1]
	
	dtheta1_dt=theta2
	dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return(dtheta_dt)

#damping coefficient
b=0.05
#gravity in m/s^2
g=9.81
#length of the pendulum in m
l=1
mass of the ball in kg
m=1
#counter variable used for naming later
ct=1

#initial condition
theta_0=[0,3]

#creating 200 time value array between 0 and 20 seconds
t=np.linspace(0,20,200)

#integrate the equation
theta=odeint(model,theta_0,t,args=(b,m,g,l))

#plotting the graphs
plt.figure()
plt.plot(t,theta[:,0],label='Angular displacement')
plt.plot(t,theta[:,1],label='Angular Velocity')
plt.legend(loc='best')
plt.xlabel('Time in seconds')
plt.show()

Now for Part 2 of the objective, I have used trigonometry to calculate the positions of the pendulum for each point of time, used the marker size command to create the sphere and finally compiled 200 images corresponding t each position of the pendulum. And using imagemagick, stitched the images together to form the animation of the pendulum. Following is the detailed code for part 2:

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


#function that returns dz/dt
def model(theta,t,b,m,g,l):
	theta1=theta[0]
	theta2=theta[1]
	
	dtheta1_dt=theta2
	dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return(dtheta_dt)

#damping coefficient
b=0.05
#gravity in m/s^2
g=9.81
#length of the pendulum in m
l=1
mass of the ball in kg
m=1
#counter variable used for naming later
ct=1

#initial condition
theta_0=[0,3]

#creating 200 time value array between 0 and 20 seconds
t=np.linspace(0,20,200)

#integrate the equation
theta=odeint(model,theta_0,t,args=(b,m,g,l))

#for creating images of the pendulum corresponding to each position
for v in theta[:,0]:
	posx=0-l*math.sin(v)
	posy=1-l*math.cos(v)
	ct=ct+1
	filename= 'output%05d.png'%ct
	plt.figure()
	plt.xlim(-2,2)
	plt.ylim(-2,2)
	plt.plot([0,posx],[1,posy])
	plt.plot([-2,2],[1,1],'-y',linewidth=10)
	plt.plot(posx,posy,'o',markersize=30)
	plt.savefig(filename)

Errors Faced:

The error was occuring as I was using the array directly inside the sin/cos function which is not possible. Also, the posx and posy are calculated wrong which I figured out soon and changed before compiling.

Solution: Used a for loop instead where the local variable contained each value of `theta` in int format for every iteration

Results and Discussion:

From the graph, it is clear that both the angular displacement and velocity follows the oscillatory trajectory which can also be seen on the actual animation of the pendulum as shown below.

The animation for the pendulum can be found below: