Week 3 - Solving second order ODEs

Second Order ODE -Simple Pendulum

Objective:

To write a program that solves the given ODE. This Ode represents the equation of motion of a simple pendulum with damping .The program also should contain the transient behaviour of a simple pendulum and to create an animation of its motion.

Solution Procedure:

ODE is used to describe the transient behavior of a system. Eg. Simple pendulum.The pendulum movement depends on newtons second law .We get a ODE which describes a position of the ball with respect to time.

Where g =gravity in m/s^2

L =length of the pendulum in m

m = mass of the ball in kg

b = damping coefficient

Given L =1 meter

m =1kg

b =0.02

g =9.81 m/s^2

Steps involved in coding:

1. Importing the different libraries such as numpy ,from scipy library the odeint inbuilt function , matplotlib for plotting the graphs and math library.
2. Constructing a function called model in which two initial conditions are given such as

dtheta1_dt = theta2

dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1)

3. Giving the constants such as length(l) , mass(m), acceleration of gravity(g) and damping coefficient.
4. We are solving two ODE, so we need to give the initial conditions for both of them.

theta_0 = [displacement, velcoity]

5. Next step is to insert the time array using linspace command.
6. Solving the ODE using odeint from scipy library where in this it will integrate the ODE over the timestep.

Here we need to give the arguments such as model which is defined as a function for ODE system where we return the dtheta_dt.

Also, theta is a solution vector, where 1^st column displays the displacement and the second column displays the velocity values.

7. .Finally, we need to plot the values for displacement and velocity over the time(t) and given the legends for representing in that graph.

   import numpy as np
   from scipy.integrate import  odeint
   import matplotlib.pyplot as plt
   import math
   
   
   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.02
   g =9.81
   l=1
   m=0.1
   
   #Initial condition
   theta_0 = [0,5]
   #time_points
   t=np.linspace(0,10,150) 
   theta =odeint(model,theta_0,t,args = (b,g,l,m))      #odeint function will integrate the ODE or  systems of ODE
   # model is a function that represents the ode system
   
   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}\theta_2-\frac{g}{l}sin\theta_1$')
   plt.xlabel('time')
   plt.ylabel('Plot')
   plt.legend(loc ='best')  
   plt.show()
   
   
   X= theta[:,0]
   ct=1
   for T in X:
   	x1=0
   	y1=0
   	x2=l*math.sin(T)
   	y2=(-l*math.cos(T))
   	file = 'test%05d.png'%ct
   	plt.figure()
   	plt.plot([-1.5,1.5],[0,0],"black",linewidth=5)
   	plt.plot([x1,x2],[y1,y2],"red",linewidth=4)
   	plt.plot(x2,y2,'o',markersize =12)
   	plt.xlim([-1.5,1.5])
   	plt.ylim([-1.5,1.5])
   	plt.savefig(file)
   	ct=ct+1
   ​

   Error-1 : dheta_dt is not defined previously. So its showing the error.

   Error 2: The constants are not declared in the arguments so its showing the error. To fix this, I had added the constants like b,g,l,m in the arguments.

   Results and Discussion:

   Hence the output gif file is uploaded to you tube.

   You Tube Link: (for gif file)

   https://youtu.be/IrhHOLk2yfQ

   Simple pendulum animation(not gif file)

   References:

   1. https://skill-lync.freshdesk.com/support/solutions/articles/43000511246-how-to-stitch-images-in-sequence-to-get-smooth-animation-
   2. http://people.uncw.edu/hermanr/mat361/ODEBook/SecondOrder.pd