Week 3 - Solving second order ODEs

Aim –

Solving a second order Ordinary Differential Equation (ODE) using Python

Objective –

1. To write a program which can solve a second order Ordinary Differential Equation (ODE).
2. To plot the displacement and velocity graph for a simple pendulum.
3. To create an animation for a simple pendulum motion.

Problem statement –

The second order differential equation used for a simple pendulum is as following -

`(d^2 θ)/(dt^2 )+ b/m (dθ)/dt+g/L sin⁡θ=0` (i)

Where,

b = damping coefficient = 0.05 Ns/m;

m = mass of bob = 0.1 kg;

g = gravitation acceleration = 9.81 m/s²;

L = length of the string = 1 m;

Above differential equation can be solved if we convert it into a first order differential equation.

Let, `θ = θ_1` (ii)

and, `(dθ_1)/dt = θ_2` (ііi)

then, `(d^2 θ_1)/(dt^2 )= d/dt( (dθ_1)/dt)= (dθ_2)/dt` (iv)

so the eq. (i) becomes,

`(dθ_2)/dt = - b/m θ_2 - g/L sinθ_1` (v)

hence eq. (iii) and eq. (v) can be written as a matrix of `(dθ)/dt`.

`d/dt [θ_1,θ_2]=[θ_2,-b/m θ_2-g/L sinθ_1]`

Problem Statement

Write a program in Python that will simulate the pendulum motion. Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0..

use,

• L=1 metre,
• m=1 kg,
• b=0.05.
• g=9.81 m/s2.

Program Outline

1. Importing required libraries matplotlib.pyplot, math, numpy, scipy.integrate .
2. Creating a function named model that takes 6 imput arguments (theta,t,b,g,l,m) and returns a list theta containing the values of displacement and velocity.
3. Defining the inputs as provided in the problem statement.
4. Setting up the initial condition provided in the problem statement (when theta is zero, velocity is 3 rad/sec).
5. Creating a time points list using np.linspace command.
6. Then the odeint function is called which takes our predefined model function, initial conditions and the input arguments and solves the equation, finally returns the list dtheta_dt which contains various values of dispacement and velocity between a period of 0 to 20 seconds.
7. In order to get the variation of displacement and velocity, a subplot is made using plt.subplot command that displays two plots of diplacement v/s time and velocity vs time and saved in a figure using plt.savefig() command.
8. For comparing both the plots, both the plots are plotted in one with necessary legents and saved as a seperate figure using plt.savefig() command.
9. In order to create the animation of pendulum for above said problem statement, for each value of displacement obtained by solving the ODE, a plot is made and saved in the working directory with specific name. for loop was used to repeat this process.
10. To create animation of these plots, ImageMagick tool was used which was activated using command promt that converted these .png files to .gif files. (The code to this has been displayed below).

Program

"""
	Solving Second Order ODE in Python
	by Saurav Sankhe
"""

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

 
# Defining function that returns theta array/list
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

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

# Initial Conditions
theta_0 = [0, 3]

# Time Points
t = np.linspace(0,20,150)

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

# Creating subplots of displacement and velocity
figure, (ax1, ax2) = plt.subplots(2)
ax1.plot(t,theta[:,0])
ax1.set_title('Displacement v/s Time')

ax2.plot(t,theta[:,1])
ax2.set_title('Velocity v/s Time')
figure.tight_layout()
plt.savefig('Displacement and Velocity plot')
plt.show()

# Creating a combined plot of Velocity and displacement
plt.plot(t,theta[:,0],'b-',label=r'$frac{dtheta_1}{dt}=theta_2$')
plt.plot(t,theta[:,1],'r--',label=r'$frac{dtheta_2}{dt}=-frac{b}{m}theta_2-frac{g}{L}sintheta_1$')
plt.legend(loc='best')
plt.savefig('Combined plot')
plt.show()

# Creating Plots for animation of Pendulum using for loop
ct = 1
for i in theta[:,0]:
	x_start = 0
	y_start = 0

	x_end = l*math.sin(i)
	y_end = -l*math.cos(i)

	filename = 'test%05d.png'% ct
	ct = ct + 1

	plt.figure()
	# Plotting base of pendulum
	plt.plot([-1, 1],[0, 0],linewidth=3,color='black')
	# Plotting string
	plt.plot([x_start,x_end],[y_start,y_end],linewidth = 3,color='blue')
	# Plotting bob
	plt.plot(x_end,y_end,marker='o',markersize = 20, color='red')
	plt.ylim([-1.25,0.25])
	plt.savefig(filename)

# End of the Program

Command Prompt code to create animation

Output Plots

From the above graph, it is clear that the initial condition, which was specified (at zero angular dispacement, angular velocity is 3 rad/sec) has been satisfied. This gives us a clear idea that the plot generated is correct.

Combined plots of Velocity and Dispacement

Animation link :

https://drive.google.com/file/d/1PTUQp8Lk45VmYKjY4EiW0IDd1HIA3Re4/view?usp=sharing

CONCLUSION

As per the initial conditions provided, the graph for the displacement starts from 0 whereas that for the velocity starts from 3m/s.

The velocity of pendulum is maximum at the mean position and 0 at the amplitude.

With every oscillation, the amplitude of oscillation decreases. For higher damping coefficient, the amplitude would decrease at a faster rate and the pendulum would come to rest early.