Week 3 - Solving second order ODEs

OBJECTIVE: To write a program that solves ODE which represents the equation of motion of a simple pendulum with damping and simulates the pendulum motion.

PROBLEM STATEMENT:  To simulate the pendulum motion if  damping coefficient,  mass of the bob,  length of the string and acceleration due to gravity  and the time for which motion continues is known . Here, It is given that :

 length of the string, L = 1m, 

 acceleration due to gravity, g=9.81 m/s^2,  

 damping coefficient,  b =  0.05

 mass of the bob, m = 1 kg

The motion is studied for the first 20 seconds and the angular position of the pendulum at t = 0 is zero degree and the angular velocity is given to be 3 radian/second.

SOLUTION PROCEDURE :

# importing libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint         # library for solving differential equations

b = 0.05                                   # Damping coefficient 
g= 9.81                                    # acceleration due to gravity
l = 1                                      # length of pendulum string in m
m = 1                                    # mass of the bob in kg

# function to determine theta array
def model(theta,t,b,g,l,m):
    theta1 = theta[0]
    theta2 = theta[1]
    dtheta1_dt = theta2                                   # 1st ODE calculation
    dtheta2_dt = (-b/m)*theta2 - (g/l)*math.sin(theta1)   # 2nd ODE calculation
    dtheta_dt = [dtheta1_dt, dtheta2_dt]                  # theta array 
    return dtheta_dt

#initial condition
theta_0 = [0,3]
# time array
t = np.linspace(0,20,100)

# odeint is a function to integrate differential equations
theta = odeint(model, theta_0, t, args = (b,g,l,m))

# plotting commands
plt.plot(t,theta[:,0],"b-",label = r'$\frac{d\theta1}{dt} = \theta2$')
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.show()

# code for saving images for aniamation 
x0 = 0
y0 = 0
ct = 0
for i in theta[:,0]:
    x = l*math.sin(i)
    y = l-l*math.cos(i)
    filename = "p" + str(ct) + ".png"
    ct = ct+1
    plt.figure()
    plt.plot([-0.75,0.75],[2,2],"black", linewidth = 5)
    plt.plot([x0,x],[2,y])
    plt.plot(x,y, marker = 'o', markersize = 20)
    plt.xlim(-2,2)
    plt.ylim(-2,3)
    plt.savefig(filename)

For solving ODE and making graph--

1. Libraries used for solving ODE such as scipy, numpy are imported along with the libraries  math  and matplotlib library and  imported the odeint()  function for integrating the ODEs.
2. Then defined the values of b,g m,L and the initial conditions that were given at t = 0.
3. The line 24 of the code  divides the time interval from 0 to 20 seconds [0,20]  into 100 parts using the function linspace() which is included under numpy module.
4. After that a function odeint() is  called  which takes some input parameters and performs internal calculations to solve the differential equation. The value that it returns is stored into a variable named  theta. The variable theta is a list which contains 100 lists in it in which every list has two elements that is angular position and angular velocity as a function of time. 
5. The parameter model that is passed into odeint() is a user defined function which tells about what the system is. This function is used for defining the two ODEs and it returns an array whose  values have to be integrated. 
                                 
   
6. After that plotting commands are used to plot the angular position and angular velocity of the pendulum with respect to time to obtain the graph.

For making the animation--

1. One end of the string which is connected to the base is fixed and is considered as origin.  Therefore the coordinates of this point are (0,0).
2. The coordinates of the other end of the string are found out using the following equation using the for loop which takes the values of  angular position from theta array.                                             x = L* sin(theta)             y = L - L* cos(theta)
3. A variable ct is used to store the plots in sequence i.e. 0,1,2,3… and so on.
4. The variable filename stores the name and extension in which the plot is going to be saved.                       For eg -:  when ct =1, filename = p1.png. 
5. Then plt.figure() command is used to create a new blank plot in each iteration of the loop. 
6. The line 45,46 and 47 of the code are for making the pendulum figure and its position with time.
7. Finally, the savefig() command is used to save the figures with the name as “filename” in a specified location.

ERRORS FACED : No major errors faced.

RESULTS AND CONCLUSIONS:

This program successfully plots the angular position vs time and angular velocity vs time graph. Since the starting point of the angular displacement plot is 0 and that of angular velocity plot is 3 which is the same as the given initial conditions, this means that the plot is correct.

From the animation, it is concluded that the oscillations die down after some time  and the amplitude of the oscillations continuously decreases with time due to damping.

Link for aniamation -

https://youtu.be/4xCxQSTwzzc

GRAPH