Week 3 - Solving second order ODEs

Second order ODE for solving pendulum with Python code and derived its transient behavior.

Title: Write a program (Python and MATLAB) to simulate the transient behavior of a system as an example of pendulum.

Objective:

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the 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.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.

Your objective is to write a program that solves the following ODE. This 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.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

Use,

L=1 meter,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Python code:

#Solving second order ODE ,A simple example is a pendulum

# this code is writen by Md.Rubayat Alam Riad
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]
	dtheta_dt = theta2
	dtheta2_dt = -(b/m)*theta2 - (g/1)*math.sin(theta1)
	dtheta_dt = [dtheta_dt,dtheta2_dt]
	return dtheta_dt

b=0.05
g=9.81
l=1
m=0.1

# initial condition 
theta_0 = [0,5]

# time points 
t= np.linspace(0,10,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\theta2}{dt}=-\frac{b}{m}\theta2-frac{g}{L}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

Output: