Week 3 - Solving second order ODEs

AIM: Write a program in PYTHON that solves the second order ODE of pendulum and animate the motion of pendulum.

Introduction

In Engineering, ODE (ordinary differential equation) is used to describe the transient behavior of a system. A Pendulum is simplest example of the ODE.

The way the pendulum moves depends on the Newtons second law. When this law is written down foe pendulum, we get a second order ODE that describes the position of the "ball" w.r.t time which as follow:

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

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.

PYTHON have odeint command to solve the ODE. But for the second order ODE first we have to convert into first order ODE.

So, convert our second order ODE in first order

First, we plot the transient behavior of pendulum than we create the animation of the motion of the simple pendulum.

Step-wise explanation of code for plot the graph of transient behavior of pendulum:

1. First, we import useful library and odeint from scipy module.
2. We have to define the function ODE_func in that the two ODE equation is written and value of some constant is given like damping coefficient(b), Length of the pendulum(L), mass of the ball(m) and gravity(g).
3. We have to assign value of theta at initial level in matrix form because it contains 2 values which are angular displacement and angular velocity.
4. We have to give the time span and convert into this span into equal interval with the use of linspace command.
5. We have to use odeint command to solve first order ODE with the help of function ODE_func with give the angular displacement and angular velocity at any instant of time.
6. Now the plot the graph of angular displacement of ball with respect to time and angular velocity with respect to time.
7. We plot these two graphs in figure so we use legend command to easy understand of angular displacement and angular velocity curve.
8. We use title command give the title of the graph.

For simulation of motion of pendulum, we have to add some steps

1. First for the show position of pendulum in different figure and for saving that we use counter.
2. We use the for loop for the show change in displacement with respect to time.
3. We assign the coordinate of string at fixed end (x₀, y₀) = (0,0)
4. We need to calculate the coordinate of free end of string (x₁, y₁) according to the value of angles which depends in time.
5. We have to plot line for pivot the string in it.
6. Then we have to plot another line for showing equilibrium position of ball.
7. Now plot the string of the graph.
8. We have to show the ball in animation so we use the marker command to plot the sphere.
9. We have to give the proper the value of X and Y axis.
10. For save every position on the ball we use name as counter command value.
11. I stich all the images in powerpoint to make animation.

PROGRAM FOR THE PLOT TRANSIENT BEHAVIOR AND ANIMATION OF SIMPLE PENDULUM

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

#defining function to give the ODE equation
def ODE_func(theta,t):

	b=0.05    #damping coeffient
	l=1       #length ofpendulm in m
	g=9.81    #gravity in m/s^2
	m=1       #mass of ball in kg
	theta1=theta[0]
	theta2=theta[1]
	dtheta1_dt=theta2
	dtheta2_dt= -(b/l)*theta2-(g/l)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return dtheta_dt

#initail condition
theta_0=[0,3]

#time span
t=np.linspace(0,20,300)

theta=odeint(ODE_func,theta_0,t)

plt.plot(theta[:,0],label='Angular Displacement')
plt.plot(theta[:,1],label='Anglar velocity')
plt.legend()
plt.title('Transient behaviour of simple pendulum')
plt.savefig("transient behaviour.png")
plt.show()

#code for the images at all the time instant
ct=1
for i in theta[:,0]:
	x0=0
	y0=0
	l=1.5 #length of string

	x1=l*math.sin(i)
	y1=-l*math.cos(i)
	
	filename = '%04d.png' %ct
	ct=ct+1

	plt.figure()
	plt.plot([-1,1],[0,0],'k-',linewidth=10)
	plt.plot([0,0],[0,-2],'k--')
	plt.plot([x0,x1],[y0,y1],'b-',linewidth=4)
	plt.plot(x1,y1,'ro',markersize=30)
	plt.xlim([-1.5,1.5])
	plt.ylim([-2.5,0.5])
	plt.savefig(filename)

	

GRAPH OF TRANSIENT BEHAVIOR OF SIMPLE PENDULUM

ANIMATION OF THE SIMPLE PENDULUM