Solving second order ODEs using PYTHON

Aim:To write a Python Program for solving the given second order ODE.

This ODE represents the equation of motion of a simple pendulum with damping.

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.

Problem Statement:Write a program in Python that will simulate the pendulum motion,

given,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

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

Introduction:A simple pendulum consists of a relatively massive object hung by a string from a fixed support. It typically hangs vertically in its equilibrium position. The massive object is affectionately referred to as the pendulum bob. When the bob is displaced from equilibrium and then released, it begins its back and forth vibration about its fixed equilibrium position. The motion is regular and repeating, an example of periodic motion.

Python Program:

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

def model(theta,t,b,g,l,m):        #Function is defined for converting the second order ODE to linear ODE so that calculation willbe easier
	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=.05   #Damping Coefficient
g=9.81  #Acceleration due to gravity (in m/s^2)
l=1     #length of pendulum(in metre)
m=1     #mass of ball(in kg)

#initial condition
theta_0=[0,3]         #First value is displacement and second value is velocity

#time points
t=np.linspace(0,20,200)    #time is defined bw 0 to 20sec and there are 200 points in between

#solve ODE
theta=odeint(model,theta_0,t,args = (b,g,l,m))  #ODE function to calculate the values

#plot results
plt.plot(t,theta[:,0],'b-',label='Displacement')  #Plotting is done with respect to time first displacement is plotted and then velocity is plotted
plt.plot(t,theta[:,1],'r--',label='Velocity')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Plot')
plt.show()


#For Animation
ct=1
for theta in theta[:,0]:                      #Animation is made with the help of positions obtained in solution
	x0=0
	y0=0
	x1=l*math.sin(theta)
	y1=-l*math.cos(theta)
	filename=str(ct)+'.png'
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],"blue",linewidth=5)
	plt.plot([x0,x1],[y0,y1],"orange")
	plt.plot(x1,y1,'o',color="green",markersize=8)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)
	ct=ct+1

Procedure:1.First of all important modules are imported for the given calculations.

2.Now the values which are given are assumed and assigned.

3.After that the function is defined so that second order ODE is converted to linear differential .Basically to make calculation easier.

4.Then with the help of ODE function the values are calculated with the given initial conditions.

5.Then plotting and animation is done.

Result: The plot obtained is :

Animation of Pendulum:

From the graph it can be easily visible that the values are decreasing with respect to time and so this is damping effect.