Solving second order ODEs using Python

Objective: To write a program in python that solves the second order ODE representing the equation of motion of a simple pendulum.

Equation of Motion representing Pendulum:

Here,

g, gravity =9.81 m/s2,

L, length of the pendulum = 1 m,

m, mass of the ball = 1 kg,

B, damping coefficient 0.05.

Calculation:

To solve the ODE of order 2, it has to be converted into a matrix containing two ODE of order 1.

In equation 5, term 1 denotes angular displacement and term 2 denotes angular velocity.

Initial Conditions:

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.

Program:

''' Simple pendulum-ODEs'''

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

#function to define the ODE
def model(theta,t,b,g,l,m):
    theta1=theta[0] 
    theta2=theta[1]
    dtheta1_dt=theta2 #displacement
    dtheta2_dt=(b/m)*theta2-(g/l)*math.sin(theta1) #velocity
    dtheta_dt=[dtheta1_dt,dtheta2_dt]
    return dtheta_dt

#inputs    
b=0.05 #damping factor
g=9.81 #gravity
l=1 #length of the pendulum
m=1 #mass of the ball

theta_0=[0,math.pi/2] #angle between the initial and first position


t=np.linspace(0,20,150) #time vector

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

#plotting
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.xlabel('plot')
plt.ylabel('time')
plt.legend(loc='best')
plt.show()

#animation
ct = 1
for i in (theta[:,1]):
    x0=0
    y0=0
    x1=l*math.sin(i)
    y1=-l*math.cos(i)
    
    filename=str(ct)+'.png'
    ct=ct+1
    
    print(filename)
    plt.figure()
    plt.plot([-1,1],[0,0],color="red",linewidth=4)
    plt.plot([x0,x1],[y0,y1])
    plt.plot(x1,y1,'o',color="green",markersize=20)
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,1.5])
    plt.savefig(filename)

Results:

Above plot denotes the relation between the X-axis Angular displacement and Angular velocity against Y-axis time.

The above image shows the motion of the pendulum. Complete animation can be viewed at https://youtu.be/m3fgZ6jsOss