Week 3 - Solving second order ODEs

Objective: To solve ODR of a simple pendulum using python.

Problem: It is necessary to solve the ODE for the given problem. It is a manual process of solving ODR and time consuming rendering. It will be very easy if it can be written in computer understandable code.

Theory: A function was defined to calculate the angular deflection and velocity of a simple pendulum at different time intervals. A second-order differential equation defining the motion of a pendulum

b = damping.

m = mass.

g = gravity.

l = length.

Problem: It is necessary to solve the ODE for the given problem. It is a manual process of solving ODR and time consuming rendering. It will be very easy if it can be written in computer understandable code.

Theory: A function was defined to calculate the angular deflection and velocity of a simple pendulum at different time intervals. A second-order differential equation defining the motion of a pendulum and defined in the function.

Method:

Import math, matplotlib, numpy, image, glob libraries
Defining 1st order ODR.
Passing initial values.
b = 0.05
g = 9.81
l = 1
m = 1
Satisfying the initial condition for theta1 and theta2.
Transition (t) for 20seconds and divide by 50 times.
Solving ODEs using odeint calls defined ODEs.
Plotting angular displacement and velocity vs time.
Graphing the state of a simple pendulum and saving the images.
Stitch all the images together to create a GIF of a simple pendulum motion.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

python code:-

'''
Simulation of motion of pendulum 
by shiek hafeez 
'''
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#function that return dz\\dt
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

#Input
b=0.02  # damping coefficient.
g=9.81  # gravity in m/s2
l=1     # length of the pendulum in m
m=1   # mass of the ball in kg

# Initial condition
theta_0 = [0,3]

# Time points
t = np.linspace(0,20,150)

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

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


# Animation
THETA3=[]
for theta3 in theta[:,0]:
	THETA3.append(theta3)

ct=0
for THETA4 in THETA3:
	x0=0
	y0=0

	x1=l*math.sin(THETA4)
	y1=-l*math.cos(THETA4)

	filename='Pendulum%05d.png'%ct
	ct=ct+1

	plt.plot([-0.2,0.2],[0,0])  # Horizontal support or pendulum
	plt.plot([x0,x1],[y0,y1])   # String of pendulum
	plt.plot(x1,y1,'-o')        # Bob tied at the end 
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1])
	plt.title('Motion of Pendulum')
	plt.savefig(filename)
	

Explanation of code:

1.To solve ode, imort odeint module from scipy module.

2.Break down the second order ODE into two first order ODE,the first ode will give the values of displacement and second ODE will give the values of velocity.

3.Define the initial condition to solve ODE as \'theta_0 = [0,3]\'.

4.To obtain smooth motion,divide time point between 0 to 20 sec using linspace command.

5.TO solve ODE use \'ODEINT\' and save values in variable \'theta\'

\'theta = odeint(model,theta_0,t,args =(b,g,l,m))\'

6.First coloumn of matrix represents the values of displacement hence extract values of displacement and add it in null matrix THETA3 using append command.

7.To show the motion of pendulum for every value of THETA3 extract values of THETA4 using for loop.

Errors:-

1.The error was in the library module i solved this by correcting scipy module

2. The error was i forgot to use math library while using mathematical function sin. solved this by using math.sin

Results: The graph below shows the variation of Angular displacement and velocity of pendulum w.r.t time.

Blue line—Displacement curve.

Redline—Velocity curve.

Below this is gif animation

youtube link:-https://youtu.be/OqPUbESldas

Conclusion:Hence ODE of simple pendulum is solved using python and respective plots were plotted with respect to time.