Week 3 - Solving second order ODEs

Objective

To simulate the pendulum motion for 20seconds.

Description

Ordinary Differential Equation (ODE) is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on Newton's 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.

Governing Equation

`(d^2 θ)/dt^2 +(b/m)(dθ)/dt +(g sin⁡θ)/l=0`

Here,

g, gravity (in m/s2) = 9.81

l, length of the pendulum (in m) = 1

m, mass of the ball (in kg)= 1

b=damping coefficient= 0.05

Code

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

# Below mentioned function is for returning 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
# given values
l=1 # length of string in meter
m=1 # mass in kilogram
b=0.05 # damping coefficient
g=9.81 #acceleration due to gravity in m/s^2


# initial conditions for solving ordinary differential equation
theta_0 = [0, 3]

# time points array
t=np.linspace(0,20,300)

# ODE solution

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

# ploting the results

plt.plot(t, theta[:,0], 'b--', label=r'$\frac{dtheta_1}{dt}=\theta_2$')
plt.plot(t, theta[:,1], 'r--', label=r'$\frac{dtheta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')

plt.ylabel('plot')
plt.xlabel('time (in second)')
plt.legend(loc='best')
plt.show()

# creating animation by using Artist Animation

fig = plt.figure(figsize=(5,5), facecolor='w')
ax=fig.add_subplot (1,1,1)

lns=[]

for i in range (len(theta)):
	#string is being plotted in below mentioned command
	ln, = ax.plot([0, np.sin(theta[i,0])],[0, -np.cos(theta[i, 0])], color='k', lw=2)

	#bob is being drawn below
	bob, = ax.plot(np.sin(theta[i,0]), -np.cos(theta[i,0]), 'o', markersize=20, color='b')

	# simulation is being plotted here

	tm=ax.text(-0.9, 0.25, 'Time=%.1fs' %t[i])
	lns.append([ln, bob, tm])

ax.set_aspect('equal', 'datalim')

# Below commands for saving animation
ani = animation.ArtistAnimation(fig, lns, interval=50)
plt.show()
fn='Pendulum_Animation'
ani.save(fn+'.mp4', writer='ffmpeg', fps=1000/50)

Code Description

• Line 1-5: Different libraries are imported.
• Line 9-15: Model function is defined for returning dtheta/dt. Here, order of ODE is reduced to solve it.
• Line 17-20: Various input parameters are assigned to variables.
• Line 24: Initial conditions are defined.
• Line 27: Time points are created using linspace command.
• Line 31: Function odeint is invoked for solving the ODE. Function takes various parameters such as model, initial conditions, time points and arguments.
• Line 35-41: Here results (angular displacement and angular velocity) are plotted using plot command.
• Line 45-68: Here animation is created using ArtistAnimation command.

Plot

Animation

For animation, refer to https://youtu.be/kcI9I4nM_Tk.

Conclusion

With the help of this code, angular displacement and velocity of the pendulum can be determined at any time point within a time frame of 20 seconds. This duration can be extended by making minor changes in code. Time duration after which pendulum will come completely to rest can also be determined.