Solving ODEs - Pendulum

The objective of this project is to study the behaviour of a transient system. i.e. A simple pendulum, by solving the second order ODE that is used to define the system.

The way the pendulum moves is defined by the second law and the equation of motion can be written as:

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

The values of these are to be assumed as:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

The second-order ODE can be written as:

`(d^2 θ)/(〖dt〗^(2 )) +b/m dθ/dt+g/L sin⁡θ=0`

Let us assume the displacement of the pendulum to be θ1 and hence:

θ = θ1 and `(dθ1)/dt=θ2`

`θ2` can be assumed to be the velocity of the pendulum.

By using the above relations, we can write the equation as:

`(dθ2)/dt= -b/m θ2-g/l sin⁡θ2`

Now we have two equations with two variables:

`d[■(θ1@θ2)]=d[■(θ2@-b/m θ2-g/l sin⁡θ2 )]`

By solving the ODE, we can get the values for the displacement and velocity of the pendulum using which we can plot the motion.

The code written for solving this equation can be seen as:

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

def model(theta,t,b,g,l,m):
	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

b = 0.05
g=9.81
l=1
m=1

#initial condition for displacement and velocity
theta_0= [0,3]

#time points 
t= np.linspace(0,20,150)
ct=1

#solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))
i=0

while i<150:
	filename = str(ct) + '.png'
	plt.figure()
	plt.plot(theta[i,0],1,'bo',markersize = 15)
	plt.plot([theta[i,0],0],[1,2])
	plt.plot([-2,2],[2,2])
	plt.plot([0,0],[0,2],'--')
	plt.xlim([-4,4])
	plt.ylim([0,4])
	plt.savefig(filename)
	i+=1
	ct+=1

plot results
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta2-\frac{g}{L}\theta1$')

plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

After importing the plotting, SciPy, math and NumPy libraries a function for defining the ODE is written. The model function takes the time interval in which the pendulum motion has to be seen and the constants (b,g,l,m) as inputs. Inside the function, the second order ODE which is used to define the pendulum motion is simplified as seen above. And an array containing the displacement and velocity is returned.

The initial conditions and the constants have been defined as seen in the code. The time interval in which the pendulum motion has to captured is set as 0 to 20 seconds. After the time interval has been set, the differential equation can be solved using the ODEINT function. The model function, initial conditions, time interval and arguments have been given as inputs to the ODEINT function.

After the equation has been solved, we can plot the displacement. At each timestep the displacement has been isolated from the solution array and it can be plotted. Since the displacement is a point a spherical marker is used to plot it. The string of the pendulum and the zero point can be plotted as lines.

After plotting the displacement, we can then proceed to plot the displacement and velocity with respect to time. We can see that the motion is damped at the end of the time interval due to the damping coefficient given.

The animation attached shows the motion of the pendulum and the displacement-velocity can be plotted as below.