Week 3 - Solving second order ODEs

AIM :

Develop a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping and simulates an animation for the simple pendulum.

Principle: -

The pendulum moves depending 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. A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with the periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained.

Solution procedure: -

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.

we used,

L=1 meter,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulation of the motion is done between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.

CODE -

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

# 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

b = 0.02
g = 9.81
l = 1.0
m = 0.1
ct = 1
# intial condition
theta_0 = [0,3]

# Time plots
t = np.linspace(0,10,150)

#Solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))
angles = theta[:,0]

# 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}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('Time')
plt.legend()
plt.show()

# Plotting 
for angle in angles:
    x1 = 5
    y1 = 5
    x2 = x1 + l * math.sin(angle)
    y2 = y1 + l * math.cos(angle)
    filename = str(ct) + '.png'
    ct = ct + 1
    fig = plt.figure(2)
    ax = plt.plot([x1,x2],[y1,y2])
    plt.plot(x2,y2,marker = 'o',markersize =10)
    plt.title('Simple Pendulum')
    plt.xlim([6.5,3.5])
    plt.ylim([6.5,3.5])
    plt.savefig(filename + '.png')
    plt.show()

• #defing function returns dz/dt - comments to make a program user friendly
• def model(theta,t,b,g,l,m): - function is used to make the code easy to implement,it is a group of statements that together perform a task
• theta1 = theta[0]
• theta2 = theta[1]
• dtheta1_dt = theta2 - initializing according to the equation
• dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1) - calculating dheta2/dt
• dtheta_dt = [dtheta1_dt,dtheta2_dt]
• return dtheta_dt - to return a value to the calling function
• import numpy as np - importing the numpy library to solve numbers in PYTHON.
• from scipy.integrate import odeint - used to solve Differential equations in Python
• import matplotlib.pyplot as plt - importing math plot library for plotting the graph
• import math - importing math library to do math operation
• b=0.02 - initializing damping coefficient
• g=9.81 - initializing gravity as 9.81m/s
• l=1 - initializing length of the pendulum as 1 m
• m=1 - initializing mass of the ball as1 kg

• # initial condition - comments to make a program user friendly

  • theta_0 = [0,3] - initializing the column of theta1 and theta2 as 0 and 3 angular displacement=0,angular velocity=3 rad/sec

  #time points - comments

  • t=np.linspace(0,10,150) - for simulation of the motion between 0-10 sec, linspace is used to generate 150values between 0 and 10

  #solve ODE - comments

  • theta = odeint(model,theta_0,t,args= (b,g,l,m)) - odeint is a command in scipy module it is a function to integrate the ordinary differential equation

  #plot results - comments

  • plt.plot(t,theta[:,0],'b-',label = r'$frac{dtheta_1}{dt}=theta2$') - command to plot the legend
  • plt.plot(t,theta[:,1],'r--',label = r'$frac{dtheta2}{dt}=-frac{g}{L}sintheta_1$') - command to plot the legend
  • plt.ylabel('Plot') - to label the y axis
  • plt.xlabel('Time') - to label the x axis
  • plt.legend(loc ='best')
  • plt.show() - to show the plot

  #Animation for simple pendulum - comments

  • ct=1 - to initialize ct as 1
  • for angle in theta[:,0]: - for executing the for loop till the displacement value in the array
  • x0=0 - initializing x0 as 0
  • y0=0 - initializing y0 as 0
  • x1=(l*math.sin(i)) - initializing x1 coordinate value
  • y1=-(l*math.cos(i)) - initializing x2 coordinate value
  • filename = give a filename
  • ct -giving the filename in sequence order
  • ct=ct+1 - incrementing the ct value to give more frames

  #Plotting

  • plt.figure()

  #Plotting the baseline of the pendulum

  • plt.plot(x1,y1,x2,y2'o',marker size =10) - plotting the bob of the pendulum
  • plt.xlim([6.5,3.5]) - giving the x limit
  • plt.ylim([6.5,3.5]) - giving the y limit
  • plt.title('Animation of Simple pendulum') - giving the title to the frames generated
  • plt.savefig(filename) - save the frames generated as filename

Output:

Animation Link (Google drive):

https://drive.google.com/drive/folders/1Sj_EEN0CyeLhYbmEFrtyKZOcD5i4lmwS?usp=sharing