Solving second order ODE's : Simple Pendulum in Python

Solving Second Order ODE's in Python : Simple Pendulum

1. AIM:

To solve a second order ordinary differential equation in Python.

2. GOVERNING EQUATIONS:

Fig 1: Trigonometry relations for a simple pendulum

Variables:

- θ\thetaθ is the angle made by rod with vertical

- B is damping coefficient

- m and L are mass of the bob and length of the rod / string

1. `(d^2θ)/(dt^2)+b/m(dθ)/dt+g/Lsinθ=0` Second order ODE. -------(1)

To solve this, we need to reduce it into lower order ODE by substitution. thus,

let `θ = θ_1`

So, we can say that `(dθ)/dt=(dθ_1)/dt=θ_2` (say) ------(2)

Therefore, `(d^2θ)/dt^2=(d^2θ_1)/dt^2=d/dt((dθ_1)/dt)=d/dt(θ_2)` ------(3)

Putting (2) and (3) in (1), we get :

`(d^2θ_1)/dt^2+b/m(dθ_1)/dt+g/Lsin(θ_1)=0` that is

`(dθ_2)/dt+b/mθ_2+g/Lsinθ_1=0`

Since it is reduced to lower order, by rearranging terms, we can say that, `(dθ_2)/dt=−b/mθ_2−g/Lsinθ_1`

3. OBJECTIVES OF THE PROJECT:

- To solve the second order ODE in Python with graphs for velocity and displacement.

- To simulate the motion of simple pendulum. (animation)

4. FLOW OF THE PROJECT:

• - Lower the order of the ODE using substitution method. Prepare a function for the same for code neatness purpose.
• - Set the geometrical and physics parameters of the pendulum like length, mass, damping coefficient as well as gravity.
• - Set the initial conditions and timespan for which you want to run the simulation.
• - Use the inbuilt ODE solver named 'odeint' to solve the ODE which is reduced by calling our function that is made in step 1 and obtain values of displacement and velocity and store it in an array.
• - Plot the graph of displacement v/s time and velocity v/s time.
• - Create simulation of the positions of the pendulum bob and arm using stitching of images by saving images first

5. PYTHON CODE:

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

#Function that returns dz/dt
def model(theta,t,b,g,l,m): 

	theta1 = theta[0]                                           #Refer governing equation section to understand these steps
	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                                                        #Damping Coefficient
g = 9.81	
l = 1															#Length of string 
m = 1															#Mass of the bob

#initial conditions
theta_0 = [0,3]													#Initial condition where angle displacement = 0 and velocity = 3m/s

#line plots
t = np.linspace(0,20,150)										#Project will for 2 seconds and t array will have 150 numbers equally spaced from 0 to 20

#Solving ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))                #Solving the ode from function by odeint function.
 
#plot results
plt.figure()                                                    #Open a new figure 

#Plotting time and angle displacements
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')

#plotting time and velocity
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()


val = theta[:,0]                                #Assingning all rows of column1  to val
ct = 1                                          #A counter to help with file names
for i in range(0,len(val)):
	plt.figure()                                #Open new figure with every iteration
	x1 = 0.5 - l*np.sin(val[i])                 #From trigonometry
	y1 = 0.5 - l*np.cos(val[i])
	plt.plot([0.5,x1],[0.5,y1])                 #Plotting the string of pendulum 
	plt.plot(x1,y1,'o',markersize=10)           #Plotting the bob
	plt.plot([-1,2],[0.5,0.5],linewidth = 4)    #Plotting the top support
	plt.xlim([-1,2])                            #Setting X and Y limits
	plt.ylim([-1,2])
	filename = str(ct) + '.png'
	plt.savefig(filename)
	ct = ct + 1
	

Code explaination:

• 1. Set physics parameters of the pendulum like mass, length, gravity and damping coefficient.
• 2. Set the initial conditions and timespan for which you want to run the simulation.
• 3. To use the ode solver, use odeint sovler from scipy module with main inputs as (t and theta) while you use model function (the one we created earlier) with inputs as (t,theta,l,b,g,m) and t and theta_0 as the initial conditions.
• 4. These results are stored in an array of theta
• 5. Since results contain 2 components, angle displacement and velocity, to plot both the graphs, use the (:,0) and (:,1) property of Python and plot those graphs.
• 6. To create animation, Use a for loop so that we can get X and Y co-ordinates using trigonometry, Here, 0.5 and 0.5 are the hinge points. Since they are (0.5,0.5) we have subtracted the co-ordinates of the bob from them so that proper rod of the pendulum can be seen.
• 7. Use the ct counter to make up for the file name using str(ct) and a '.png' file extension.
• 8. Use a image stitcher to stich all the images to create an animation.

6. ERRORS FACED:

- Fixed by changing val=theta[:,1] to val=theta[:,0]. Since indexing in python starts from 0.

- Also, the combining of plots is fixed my putting plt.figure() inside of the for loop instead of outside.

7. OBSERVATIONS:

- We can see that the plot of Displacement(blue color) and Velocity(red) is plotted against time. We can see the reduction of velocity and angle displacement as time increases.

8. CONCLUSION:

- The reduction in velocity and Angle displacement is due to damping coefficient introduced in the Second order ODE of a simple pendulum.

- We have successfully solved the second order ODE.

- The animation of the pendulum is here https://www.youtube.com/watch?v=DlL59E_YgEo

9. REFERENCES:

- Figure 1 - Khan Academy