Week 3 - Solving second order ODEs

AIM: To solve the 2nd order ODE and alsos to animate the motion of the pendulum with the help of Python.

Implementation:

The oscillation of the pendulum with respect to the angular displacement and the velocity is been calculated using the odeint of scipy.integrate module where the result is a form of an array that is plotted and also the position of the pendulum with respect to the time is been animated using the required code as stated below.

The ode function was developed based on the eqn

there by the equation was broke into smaller single degree eqns for simplicity as follow

theta=theta1

dtheta1/dt=theta2

d^2theta/dt=dtheta2/dt

thereby the the above equation got reduced to the form of dtheta2/dt=-(b/m)theta2-(g/l)sin(theta1)

The following if the code for the solving of the ode and also to animate the motion

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

#function that returns 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 parameters
b=0.05
g=9.81
l=1
m=1

#initial condition
theta_0 = [0,3]

#time points
t= np.linspace(0,20,200)

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

#plotresults for the ode solved
plt.figure('1')
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.savefig('disp_vel_Vs_time')

#now to animate the motion of the pendulum with a bob at the end.
x0=0
y0=0
l=1
ct=1 #counter variable 
for i in range(0,len(t)):
	x1=-l*math.sin(theta[i,0])
	y1=-l*math.cos(theta[i,0])
	
	filename=str(ct)+'.png'
	ct=ct+1
	plt.figure()
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker="o",markersize=20)
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.savefig(filename)

The above code consists of 3 parts in short.

1. ODE function to solve the above given equation.
   Here the equation is broke into individual single degree equations and finally clubbed in the form of a matrix i.e. array of angular displacement and angular velocity respectively and then returned that final value of theta.
   
   
2. Using the defined ode function to solve the the ODE and also to plot the results so obtained
   Here the input parameters are taken as damping coefficeint (b)=0.05, g=9.81m/`s^2`, mass (m)= 1kg and length of the pendulum (L)=1m, angular displacement=0,angular velocity=3 rad/sec and a time period of 20 seconds.
   Taking these inputs, the final result of array theta was found and then plotted accordingl using the matplotlib.pyplot module.
   
   
3. Now animating the motion of the pendulum.
   In order to plot the position of the pendulum it was needed to find the positions of the pendulum at each position where one position was fixed as x0=0 and y0=0. Now a for loop is initiated throughout the length of the time interval and by finding the end positions x1,y1 with respect to the angular displacement found by solving the ODE.
   Counter variable was used to count the iterations and also to help in naming the individual plots to save and later to form an animation of it.
   limits for the x and y axis were given respectively and also the marker was plotted to visualize the weight hanging in motion.

Outputs.

Plot showing the displacement/velocity variation with respect to time.

Now the animation of the pendulum in motion which was done using the plots generated by python and converted to animation using imgflip.com .

Below is the screenshot showing the plots being generated by the python code.

Errors
Most of the errors were typo's which were immediately sorted out and the syntax of the ploting i.e. for marker size. This was rectified by checking the actual syntax for increasing the size of the marker.

the right syntax for the marker is plt.plot([x1,y1],marker="o",markersize=20) , this is used in the above code.

Conclusion:

The ODE is solved using the python code and also the pendulum motiion is animated as per the requirement stated.

References:

• https://kite.com/python/answers/how-to-set-point-size-in-a-matplotlib-plot-in-python