Week 3 - Solving second order ODEs

Aim:

Write a program to simulate the transient behaviour of a simple pendulum and create an animation of its motion.

Introduction:

A simple pendulum is defined to have an object that has a small mass, also known as the pendulum bob, which is suspended from a light wire or string.

In Engineering, ODE is used to describe the transient behaviour 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.

ODE represents the equation of motion of a simple pendulum with damping

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.

Initial condition:

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.

Python Program:

# Simulation of Pendulum Motion

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


# function that returns dz/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

# Input
b = 0.05 # damping coefficient
g = 9.81 # gravity in m/s^2
l = 1    # length of the pendulum in m
m = 0.1  # mass of the ball in kg

# initial condition
theta_0 = [0,3]

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

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

# Plot Results
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend( loc = 'best')
plt.title('Second Order ODE')
plt.savefig('Second Order ODE')
plt.show()

# Animation
Theta3 = []
for theta3 in theta[:,0]:
	Theta3.append(theta3)

ct = 1
for Theta4 in Theta3:
	x0 = 0
	y0 = 0

	x1 = 1*math.sin(Theta4)
	y1 = -1*math.cos(Theta4)

	filename = 'Pendulum%05d.jpg'%ct
	ct = ct+1

	plt.plot([-0.2,0.2],[0,0]) # Horizontal support or pendulum
	plt.plot([x0,x1],[y0,y1])  # String of pendulum
	plt.plot(x1,y1,'-o')	   # Bob tied at the end
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1])
	plt.title('Pendulum Motion')
	plt.savefig(filename)
	plt.show()

Explanation of the Program:

• To solve ODE, import the odeint module from the scipy module.
• Break down the second-order ODE into first-order ODE, the first ODE will give the values of displacement and the second ODE will give the values of the velocity.
• Define the initial condition to solve ODE as 'theta_0 = [0,3]'.
• To obtain smooth motion, divide the time point between 0 to 20 sec using the linspace command.
• To solve ODE use 'ODEINT' and save values in variable 'theta' 'theta = odeint (model, theta_0, t, args = (b, g, l, m))'
• The first column of the matrix represents the values of displacement hence extract values of displacement and add them to null matrix Theta3 using the append command.
• To show the motion of the pendulum for every value of Theta3 extract values of Theta4 using the loop.
• The obtained images are merged into mp4 video format in command windows as shown in the figure below:

Result:

1. Time Plot

2. Animation of Pendulum Motion

https://youtu.be/IJVG5wHAavA

Thank you...