Week 3 - Solving second order ODEs

Objective-

Write a program that solves the following ODE by using Python. Plot the graph of Differential equation and animation of the Pendulum using Python.

Problem Statement-

ODE is used to describe the transient behaviour of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Equation-

The ODE represents the equation of motion of a simple pendulum with damping

where

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Here to solve the equation the values given

b=0.05

g=9.81 m/s2

m=1kg

l=1m

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

Code-

#Program to solve the ODE 
#By Vinod Kumar
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt 
import math

#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
#constants
#damping coefficient
b=0.05 
#acceleration due to gravity 
g=9.81
#length of the pendulum
l=1
#mass of the pendulum ball
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))

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

# For Animation of the simple pendulum
ct=1
for theta in theta[:,0]:
	#coordinate of the base point
	x0=0
	y0=0
	#coordinate of the end point of pendulum
	x1=l*math.sin(theta)
	y1=-l*math.cos(theta)

	filename='test%04d.png'%ct
	ct=ct+1
	plt.figure()

	#for Horizontal base line of pendulum
	plt.plot([-0.5,0.5],[0,0],'r',linewidth=4)
	#plotting string of pendulum
	plt.plot([x0,x1],[y0,y1],'g')
    #plotting pendulum ball
	plt.plot(x1,y1,'o',markersize=20)

    #limiting x-axis
	plt.xlim([-1.5,1.5])

    #limiting y-axis
	plt.ylim([-1.5,1])

    #Title of the animation 
	plt.title('Animation of Simple Pendulum using Python')

    #saving the files
	plt.savefig(filename)

Output-

The ODE equation Graph-

The animation of the simple pendulum

Solution Procedure-

1. Import the numpy module, scipy module, matplotlib.pyplot module & math module

2.Define the system of differential equation

3. Convert the differential equation into 1st degree of the differential equation.

4. Take initial conditions and input values of constants(b,l,m,g).

5. Here the given time span is 0-20sec so we take 200 equally spaced value between 0 to 20sec.

6. Use 'odeint' command from the scipy module that integrates the ordinary differential equations.

7. Plot the graph of Angular displacement and Angular velocity with time.

8. Use limiting value of x-axis and y-axis and the legend, label and title.

for animations

9. Take value '1' for ct.

10. Take the array 'theta' to set the values of angular displacement.

11. Define the initial coordinate of pendulum and endpoints with the mathematics equation.

12. Use FOR loop to run the code for different angular displacement.

13. Plotting the Horizontal line, String and pendulum ball.

14. Use limiting value of x-axis and y-axis and the title.

15. Use Savefig(filename) to save each position of the pendulum for each angular displacement.

Results-

1. The Graph of the Ordinary Differential Equation for the Simple Pendulum successfully plotted with the help of Python.

2. The animation for the Simple Pendulum is successfully developed with the help of python