Solving second order ODEs using Python

• AIM : To simulate the transient behaviour of a simple pendulum and to create an animation of it's motion on python
• objective this poject is about the how to calculate the simple harmonic motion in simple pendulum using python and how to make the animation in this pythone code odeint function has been used to solve the problem.
• Theory

A simple pendulum consist of a metal mass held by a string to a stationary wall.The mass oscillates because of gravity and its inertia when lifted and then left to fall.This oscillation can be modeled using second order differential equation. 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.

We consider a simple pendulum that consists of a mass m suspended by a string or rod of length L (and negligible mass). The angle it makes with the vertical varies with time as a sine or cosine.

Governing equation

The ODE which represents the equation of motion of simple pendulum with damping is

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

In mathematics, an ordinary differential equation (ODE) is a differerntial equation containing one or more functions of one independent variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable

• python code :

#importing important packages
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

#inputs
b = 0.05
g = 9.81
l = 1
m = 1

#initial condition

theta_0 =[0,3]

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

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

#plot result
plt.plot(t,theta[:,0],'b-',label='angular displacement')
plt.plot(t,theta[:,1],'r--',label='angular velocity')
plt.xlabel('time')
plt.ylabel('vectors')
plt.legend()
plt.show()
plt.savefig('displacement and velocity graph.png')

#for animation

ct = 1


for theta in theta[:,0] :
	x0 =0
	y0= 0

	x1=(1*math.sin(theta))
	y1=(-1*math.cos(theta))
	filename = str(ct)+'.png'	
	ct = ct+1	
	plt.figure()	
	plt.plot([-1,1],[0,0],"black",linewidth = 10,marker="o")
	plt.plot([x0,x1],[y0,y1],"red")
	plt.plot(x1,y1,'o',markersize = 20)
	plt.xlim([-1.5,1.5])	
	plt.ylim([-1.5,1.0])
	plt.savefig(filename)

code explanation :

The libraries required for this program are imported.Especially, in order to solve ordinary differential equation, 'scipy' module of 'odeint' solver is imported.Then using def command we create aFunction named 'model' is defined to solve those converted simple ordinary differential equations.The calculated values are stored in the array 'dtheta_dt' and 'return' command is used to return those values when the function is called in the main program. then we created a code for main function first type the inputs to be given - length(l) = 1m, mass(m) = 0.1kg, damping(b) = 0.05, gravity(g) = 9.81 we are giving the intial condition as theta[0, 3] which means that the intial displacement is 0 and initial velocity as 3 rad/sec. next we are using linspace (np.linspace) command to generate values 200 values between 0 and 20 second for the time span then odeint command used to solve 2 frist order ODE equation in the function and give the result .
Here, the endpoint of the pendulum is fixed and 'for loop' is used to plot the positions of the bob with respect to time with the calcuated values in the 'results_ode' array.The pendulum bob is created using 'markers' and 'markersize' and .The 'savefig()' command is used to get each and every plot and save images as '.png' format from the array 'filename'

Here, the values of angular displacement and angular velocity is plotted with time and the variation is shown and the respective equations also labelled in the plot.

• plot/result :

plot for time and angular velocity,displacement

link for the motion of pendulum :

https://youtu.be/5cHIo8Udloc

CONCLUSION:

Thus the program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion is coded and executed successfully by uploading the video in youtube.