Solving second order ODEs and making animation for movement of pendulum using the ODEs in python

AIM -

To write a python program to solve the second order differential equation and run a animation for pendulum

OBJECTIVE -

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

SECOND ORDER DIFFERENTIAL EQUATION -

The order of a differential equation is the highest-order derivative that it involves. Thus, a second order differential equation is one in which there is a second derivative but not a third or higher derivative.

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.

CODE SAMPLE -

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

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

b=0.2
g=9.81
l=1
m=0.1

# Initial condition
theta_0 = [math.pi/2,0]

#time points
t = np.linspace(0,10,150)

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

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('time')
plt.ylabel('plot')
plt.legend(loc='best')
plt.show()


# for animation of pendulum

ct = 1
for i in (theta[:,1]):
	
	x0 = 0
	y0 = 0
	x1 = 1*math.sin(i)
	y1 = -1*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',markersize = 20)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)

EXPLANATION OF CODE -

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

First we have imported the required modules for solving the ode problem. math and matplotlib are the common modules to generate mathematical and plot functions. numpy is used to form a array or more defining an array in proper manner. scipy module is must to solve the ODE by using odeint.

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

Above def function basically defines the theta,t,b,g,l,m that are g = gravity in m/s2, L = length of the pendulum in m, m = mass of the ball in kg, b=damping coefficient, t= time interval in sec.

The required values used in code are -

b=0.2
g=9.81[m/s2]
l=1[m]
m=0.1[kg]

theta_0 = [math.pi/2,0]

t = np.linspace(0,10,150)

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

Above the theta is defined from 90 degree to 0 at time span of 0 to 10 at 150 interval using np.linspace command. Basically, np.linspace is used to define or write the given function within prescribed range i.e., here we have time span from 0 to 10 and we see that 150 is also mentioned so what does 150 do is it generates 150 values between 0 to 10. Futher, odeint is used to solve the differential equation by using the mentioned values.

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('time')
plt.ylabel('plot')
plt.legend(loc='best')
plt.show()

So, above alogrithm is basically the plot command.

The above plot shows us the variation of two equation mainly angular velocity and angular displacement with the time.

PENDULUM ANIMATION CODE -

ct = 1
for i in (theta[:,1]):

x0 = 0
y0 = 0
x1 = 1*math.sin(i)
y1 = -1*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',markersize = 20)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)

Above, code for the animation of the pendulum. The for loop is used for theta to access the values for pendulum movement and defined the values foe x0, y0, x1, y1.

The generated images were saved and futher were converted to .gif

The link for the pendulum animation is provided below.

https://youtu.be/9MbkoIFd34w

Conculsion - By using the ODEs we can solve many engineering apllication that can be used in day to day time.