Solving second order ODEs using python

Objective : To write a program in Python that solves the second order ODE corresponding to motion of a simple pendulum and to plot its angular displacement and angular velocity wrt time & to create the animation of the motion of pendulum by simulating the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.

Description- A simple pendulum consists of a point mass m hanging on a string of length L from one support. One pulls the mass back to some starting angle θ0">θ0 and releases it than the goal is to find the angular position as a function of time. The way the pendulum moves depends on the Newton's Second Law of Motion. When this law is written down, we get a second order differential eqauation that describes the position of the ball w.r.t time. The simple pendulum is described by the below mentioned ODE-

where,

g = acceleration due to gravity (m/s^2)

L = length of the string

m = mass of the ball

b = damping coefficient

Given that

g = 9.81 (m/s^2)

L = 1 meter

m = 1 kg

b = 0.05

For angular Displacement = 0

Angular velocity = 3 rad/sec

at time t=0.

Simulate the motion between 0-20 sec

PYTHON CODE FOR SOLVING ODES AND ANIMATION

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
def ode_function(theta,t,b,g,l,m):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = - (b/m)*theta2-(g/1)*math.sin(theta1)
	dtheta_dt = [dtheta1_dt,dtheta2_dt]
	return dtheta_dt

# Inputs
# damping coefficient
b = 0.05

# acceleration due to gravity(m/s^2)
g = 9.81

# length of string (m)
l = 1

# mass of bob (kg)
m = 1

# Initial Conditions
theta_0 = [0,3]

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

# Solve 0DE
theta = odeint(ode_function, theta_0,t_span,args = (b,g,l,m))

# Plotting the graph
plt.plot(t_span, theta[:,0], 'b-',label=r'$Displacement$')
plt.plot(t_span, theta[:,1],'r--',label=r'$Velocity$')
plt.xlabel('Time (second)')
plt.ylabel('Displacement or Velocity')
plt.title('Damped Harmonic Motion of a Simple Pendulum')
plt.legend(loc='best')
plt.grid()
plt.savefig('')
plt.show()


# Iteration using for loop to create animation 
loop_counter = 1
tmp = theta[:,0]
for i in tmp: 
	x_start = 0
	y_start = 0
	x_end = l*math.sin(i)
	y_end = -l*math.cos(i)

	# Plotting for animatio
	plt.figure()
	
	# Plotting String of pendulum
	plt.plot([x_start,x_end], [y_start,y_end],linewidth=3)
	
	# Plotting Bob of pendulum
	plt.plot(x_end, y_end, 'o',markersize=20, color='k')
	
	# Plotting the rigid support through which string of pendulum is attached
	plt.plot([-0.5,0.5],[0,0],linewidth=8,color='y')
	
	# Plotting hatched vertical line for refernce
	plt.plot([0,0], [0,-1.5],'--', color='r')
	plt.xlim([-1.25,1.25])
	plt.ylim([-2, 1])
	plt.title(' Animation of transient behaviour of a Simple Pendulum')
	filename = 'test%05d.png' % loop_counter
	plt.savefig(filename)
	loop_counter = loop_counter + 1

Explanation of code -

• At first math, matplotlib, numpy and scipy modules were imported for their various functions.
• Initial inputs such as b,g,l,m are provided according to data in the question.
• An array of initial theta was created to define initial angular displacement and initial angular velocity.
• An array of time span was created to store 200 equally spaced intervals between time 0 to 20 seconds.
• To solve the second order ODE, a function named ode_function which stores the above matrix data in it.
• The given second order ODE was solved using 'odeint' command followed by proper syntax of the calling function and the result obtained was used to plot the graph of variation of angular displacement and angular velocity wrt time.
• Now for creating the animation of obtained harmonic motion of pendulum, for loop was used and a loop counter variables was initiated,the values of angular position is stored in a separate variable called 'tmp' and is iterated in for loop.
• In this execution the one end of the string of pendulum was fixed and the other end moves according to the solution of the ODE solved and was plotted for each time interval of the time span.

RESULTS :

The following plot was generated for the given program –

Conclusion: -

• The above graph shows the angular displacement and angular velocity and it changes over the time due to the effect of the damping, the values of reduce angular displacement and angular velocity with time as is clearly visible in above plot.
• The google drive link of the animation is https://drive.google.com/file/d/1Cnx_Ht0IXu5Xhc8Vzj755y43Wa8Q7-Zk/view?usp=sharing