Solving second order ODE using PYTHON

Aim :- To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

`frac{d^2theta}{dt^2} + frac{b}{m} frac{d theta}{dt} + frac{g}{L} sintheta = 0`

g = gravity in m/s2,L = length of the pendulum in m,m = mass of the ball in kg,b=damping coefficient

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.

Code :-

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

#Defining Function that returns 'dtheta_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


#Inputs
b = 0.05
g = 9.81
l = 1.0
m = 1.0

#Initial Condition
theta_0 = [0,3]

#Time Points
t = np.linspace(0,20,150)

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

#Plotting Results 

#Plotting Angular Displacement
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta_2$')

#Plotting Angular Velocity
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')

#Limiting X-axis
plt.xlim([0,20])
#Limiting Y-axis
plt.ylim([-3,3])

plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.show()


#animation
ct = 1
for theta in theta[:,0]:
	x0 = 0
	y0 = 0                                                         # x0,y0 fixed end of pendulum
	x1 =  l*math.sin(theta)
	y1 = -l*math.cos(theta)                                        
	plt.figure()
	plt.plot([-0.25,0.25],[0,0],"black",linewidth=2)                    
	plt.plot([x0,x1],[y0,y1],"red")                               
	plt.plot(x1,y1,'o',markersize=20)                             
	plt.xlim([-1,1])
	plt.ylim([-1.5,1])
	ct = ct +1
	filename = 'test%05d.png' % ct
	plt.savefig(filename)
	

Output :-

OUTPUT :- Animation of the pendulum in you tube.

Error :- Type of error getting during this project is return pendulum to it's starting position and not not having easily smooth animation, this is easily removed by using magick *.png output.gif.