Week 3 - Solving second order ODEs

Objectives:

• Write a program in python that would solve the second order ordinary differential equation of motion of a simple pendulum with damping

Theory:

A simple pendulum is a system having a mass of small diameter, also reffered as bob and a string that has a negligible mass but is strong enough to stay taut. The linear displacement from equilibrium is equivalent to the length of the arc described the bob with radius equal to the length of the string. Also the force acting on the bob when displaced from its equilibrium position is mgsin`mgsinѲ which result in a net force acting towards the equilibrium position, reffered as restoring force.

the motion of the pendulum mass is governed by the Newton's second law of motion. When this law of motion is written down, a second order Ordinary Differential Equation that describes the position of the bob w.r.t time is derived.

The derived mathematical equation 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

Solution Procedure:

The second order ODE for simple pendulum with damping is:

In order to solve equation (1), it is converted to two simultaneous ODEs of order 1 by the below method:

Variable Definition:

b = 0.05

m = 1 kg

L = 1 m

g = 9.81 m/s^2

Initial condition:

The values of position and angular velocity at initial condition, i.e at t=0 sec is defined in an array named theta_0 as displacement =0 and angular velocity = 3 rad/sec

Time array:

A time array is created using numpy module for 500 time steps between 0 sec and 20 sec.

Function definition:

A function named pendulum_system is defined to solve the two simultaneous ODEs of the first order. The function is called for the values of the following:

• Empty array to store the solution (Theta)
• time array (t)
• damping coefficient (b)
• acceleration due to gravity (g)
• String Length (L)
• Mass of the bob (m)

Solution of the ODE:

The ODE is solved by using odeint command imported from the scipy.integrate module. It takes into account the following varialbes to return the solution array:

• pendulum_system function
• initial condition (theta_0)
• time array (t)
• argument of L, m, b, g

Results and Discussions:

The solution gives us a plot of displacement Vs time and Angular velocity Vs time. The graph would vary depending on array values in the initial condition array (theta_0)

Code:

# Solving second order ODE

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

# damping coefficient
b = 0.02
# acceleration due to gravity (m/s^2)
g = 9.81
# length of the string (m)
l = 1
# mass of the system (kg)
m = 0.1

# function returning the velocity values
def pendulum_system(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


# initial conditions
theta_0 = [0, 3]

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

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

# plotting the graphs
plt.plot(t, theta[:,0], "b-", label=r'$frac{dtheta_1}{dt}=theta2$')
plt.plot(t, theta[:,1], "r--", label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{L}sintheta_1$')
plt.xlabel("time")
plt.ylabel("plot")
plt.legend(loc='best')
plt.show()

# Animation

# Counter
ct = 1

for i in theta[:, 0]:
	x0 = 0
	y0 = 0
	x3 = -1.5
	y3 = 0
	x4 = 1.5
	y4 = 0
	x1 = l*math.sin(i)
	y1 = -l*math.cos(i)
	filename = str(ct) + '.png'
	ct = ct +1
	plt.figure()
	plt.plot([x0, x1], [y0, y1])
	plt.plot([x3,x4], [y3, y4])
	plt.plot(x1, y1, 'o', markersize = 20)
	plt.xlim(-1.5, 1.5)
	plt.ylim(-1.5, 0)
	plt.savefig(filename)