Week 3 - Solving second order ODEs

INTRODUCTION

Many physical systems can be described using ordinary differential equations (ODEs), solving the ODEs means that the solution of the ODE can be used to predict the behaviour of the physical system. There are dedicated modules in the numpy module of the scipy library in python which can be used to solve an ODE.

AIM

In this challenge, we are solving an ODE which describes the motion of a simple pendulum and an animation is created which shows the motion of the pendulum.

CODES and EXPLANATION

The equation that we are trying to solve is a second-order ODE, which is given by:

`\frac{d^2\theta}{dt^2} + \frac{b}{m}\frac{d\theta}{dt} +\frac{g}{l}sin\theta =0` which is an equation which describes a damped free vibration in mechanics.

The above ODE can be solved using odeint() function from the scipy library. The explanation of the code is given below

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

The above snippet of code imports the required libraries. The numpy library is used for managing the arrays, scipy is the scientific computing library which contains tools which help in scientific computing, matplotlib is used for creating the plots and math library is used for mathematical functions.

"""
Solving the ODE
"""
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

In the above snippet of code, a function is defined which act as one of the parameter inputs to the odeint() function which solves the ODE. Here the second-order ODE is split into two different first-order ODEs and the above function returns a matrix consisting of two first-order ODEs.

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

# initial Condition

theta_0=[math.pi/2,0]

# time points
t=np.linspace(0,20,150)

#solve ode
theta=odeint(model,theta_0,t,args=(b,g,l,m))

In the above snippet b represents the damping coefficient, g is the acceleration due to gravity l is the length of the pendulum and m is the mass of the bob of the pendulum. The initial condition is an array with two entries, first is the initial angular displacement and the second entry represents the initial angular velocity. The solution is integrated up to 20 s, which is represented by t, where the time interval is divided into 150 equal intervals using the linspace() function in the numpy library.

The odeint() function returns the angular displacement and the angular velocity in the form of a matrix, where the first column consist of angular displacement and the second column consist of angular velocity.

plt.plot(t,theta[:,0],'b-',label = r'$\frac{d\theta_1}{dt} = \theta_2$')
plt.plot(t,theta[:,1],'r--',label = r'$\frac{d\theta_2}{dt} = -\frac{b}{m}\theta_2 - -\frac{g}{l}sin\theta_1$')
plt.legend(loc ='best')
plt.ylabel('Plot')
plt.xlabel('time')
plt.show()

The above snippet of code is used to plot the displacement and velocity vs time. This theta[:,0] access only the displacement terms and this theta[:,1] only access the velocity terms. The output of the above code is shown below

The code used for creating the solution animation of the pendulum is given below:

a = theta[:,0]
x0 = 0
y0 = 1
count = 0 
for i in a:
    x1 = l*math.sin(i)
    y1 = -l*math.cos(i)
    count = count + 1
    filename = 'test%05d.png' % count
    plt.figure()
    plt.plot([-1,1],[1,1],linewidth = 5)
    plt.plot([x0,x1],[y0,y1], linewidth = 2)
    plt.plot(0,1,'ro',markersize = 5)
    plt.plot(x1,y1,'bo', markersize=10)
    plt.xlim(-1.25,1.25)
    plt.ylim(-1.25,1.25)
    plt.axis('off')
    plt.savefig(filename)
    

An array is created in the first line and all the angular displacement value is stored there. The value of all those angular displacements is used to calculate the position of the bob. Here (x0,y0) represents the string attached to the roof and (x1,y1) represents the position of the bob. The x coordinate is represented by `lsin\theta`and the y coordinate is represented by `lcos\theta`. The bob position for different values of theta is calculated and plotted. All the plots are stitched together using a python code to create a gif animation which is shown below.

The code written in python is shown below:

"""
Simple Pendulum animation

@author: Abhijith S Pai
"""
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

 
"""
Solving the ODE
"""


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.05
g=9.81
l=1
m=1

# initial Condition

theta_0=[math.pi/2,0]

# time points
t=np.linspace(0,20,150)

#solve ode
theta=odeint(model,theta_0,t,args=(b,g,l,m))

plt.plot(t,theta[:,0],'b-',label = r'$\frac{d\theta_1}{dt} = \theta_2$')
plt.plot(t,theta[:,1],'r--',label = r'$\frac{d\theta_2}{dt} = -\frac{b}{m}\theta_2 - -\frac{g}{l}sin\theta_1$')
plt.legend(loc ='best')
plt.ylabel('Plot')
plt.xlabel('time')
plt.savefig('displacement.png')
plt.show()


""""
Creating the plots for animating 
"""

a = theta[:,0]
x0 = 0
y0 = 1
count = 0 
for i in a:
    x1 = l*math.sin(i)
    y1 = -l*math.cos(i)
    count = count + 1
    filename = 'test%05d.png' % count
    plt.figure()
    plt.plot([-1,1],[1,1],linewidth = 5)
    plt.plot([x0,x1],[y0,y1], linewidth = 2)
    plt.plot(0,1,'ro',markersize = 5)
    plt.plot(x1,y1,'bo', markersize=10)
    plt.xlim(-1.25,1.25)
    plt.ylim(-1.25,1.25)
    plt.axis('off')
    plt.savefig(filename)
    

 

 



 

 

CONCLUSION

The code is written in python to solve an ODE and the solution is used to create an animation of a simple pendulum.