Week 3 - Solving second order ODEs

Aim: To write a program and plot a graph that solves the ODE represents the equation of motion of a simple pendulum with damping

Theory:

ODE:An Ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The equation which we are using here is :

In the above equation,

g = gravity in m/s^2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b= damping coefficient.

Input Variables for the Animation and Plotting of the Graph

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Code:

# -*- coding: utf-8 -*-
"""
Created on Fri Jun 19 22:40:25 2020

@author: PANKAJ
"""

#Importing the required modules
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
from PIL import Image
import matplotlib.image as image
import glob


#Defining the function to return the value of 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

#Input values for the pendulum 
    
b=0.05
g=9.81
l=1.0
m=1.0


#Initial Condition
theta_0=[0,3]

#Time Intervals

t= np.linspace(0,20,200)

#Solving the differential equation
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#Plotting  the results
#Plotting the angular displacement
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta1}{dt}=\theta2$')
#Plotting the angular velocity
plt.plot(t,theta[:,1],'r-',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta2-\frac{g}{l}sin\theta1$')
plt.legend(loc='best')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.show()


#Creating the Animation
c=1  #Creating the counter loop variables
frames = [] #The images created would be appended
imgs = glob.glob("*.png") #Using glob method to extract .png files

for x in theta[:,0]:
    x0=0
    y0=0
    x1=l*math.sin(x)
    y1=-l*math.cos(x)
    #geting an image of wall as the base of the pendulum
    im = Image.open("Pendulum wall.png")
   
 
      
    filename=str(c)+'.png'
    c=c+1
    plt.figure()
    #plotting the image file and setting the axes and location
    fig, ax = plt.subplots() 
    
    ax.imshow(im,aspect='auto',extent=(-1.0,1.0,0.0,0.5),alpha=0.5)
    
    
    #PLotting the String of the Pendulum
    plt.plot([x0,x1],[y0,y1],'r')
    #PLotting the Ball of the Pendulum
    plt.plot(x1,y1,'o',markersize=20)
    
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,1.5])
    plt.title('Animation of Pendulum using Python')
    #Saving the figures with their file names
    plt.savefig(filename)
    plt.close()
    for i in filename:
    			new_frame = Image.open(filename)
    			frames.append(new_frame)
 

		# Create the frames from the extracted files from glob

# Saving into a GIF file that loops forever
frames[0].save('Pendulum Animation.gif', format='GIF',append_images=frames[1:],save_all=True, duration=50, loop=0)
print('Done')

Code Explaination:

Step1:

Importing the required modules such as numpy, scipy.integrate, matplotlib.pyplot,math,Image,Glob, matplotlib.image.

Step2:

Defining the function to return the value of dtheta_dt which takes the variables angle, time, gravitational acceleration,length,mass and damping coefficient.

Step3:

Initializing the input values for the pendulum. Setting the velocity interval between 0 to 3 rad/s and the time span from 0 to 20 sec

Step4:

Solving the differential equation and finding the value of theta using the odeint method of the scipy module which is used for the integration of the model function that was created earlier. From this we get theta which is a solution vector. It contains two columns , solution for the first ode integration gives displacement values and the second column gives the velocity values.

Step 5:

Plotting the angular displacement from the first column of theta value and angular velocity from the second column of the theta vector.

Step6:

Creating the Animation

Creating a counter loop variable c and an empty array in which the images would be appended. Using the glob method to extract .png files in the folder. Creating a ‘for’ loop for the displacement vector of theta in which initialization of the coordinate variables are done, they are the location of the bob of the pendulum. Then loading an image of wall as the pivot of the pendulum which is inserted and plotted using the subplot method of matplotlib. Then the string of the pendulum is plotted using the distance between the pivot and the bob, similarly the bob is plotted at the end of the string using the marker, and plot is set to a limit. Then the file is saved and this is saved for 200 times with each different displacement value of theta.

The frames are created from the extracted files from glob and saved into a gif file that loops forever.

Plot:

Errors:

This error occurred while appending the image in the plot. To overcome this I directly used the Image module from matplotlib and used it to subplot the image and then set the axes using the extent variable of the imshow method.

Youtube link for the Animation:

https://youtu.be/kkmd8AtKICI