Solving second order ODEs using PYTHON

Aim: Solving Second Order Differential equations using PYTHON.

Theory: A simple pendulum is one which can be considered to be a point mass suspended from a string or rod of negligible mass.

When displaced from its equilibrium point, the restoring force which brings it back to the centre is given by:

`F_"net"=−m⋅g⋅sinθ`

For small angles θ, we can use the approximation

`sinθ ≈ θ`

in which case Newton's 2nd law takes the form

`F_"net" = ma = −mgθ`

Even in this approximate case, the solution of the equation uses calculus and differential equations. The differential equation is

`(d^2 θ)/(dt^2 )+g/L θ=0`

However, the above case is assumed that there is no damping in the system. When we introduce damping, we get a modified equation as follows:

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/L θ=0` ....................................................................... (A)

Where, b = damping factor and m = mass of the bob

Above equation is of second order, so we will try to convert the above equations into smaller differential equations which will lead us to the following matrix

`d/dt((theta_1),(theta_2))=[[theta_2],[-b/m*theta_2-g/L*sin(theta_1)]]`

It is done as follows:

Let `theta = theta_1` so differentiating above equation w.r.to t we get

`d/dttheta = d/dttheta_1 =theta_2` .............................................................................. (1)

differentiating above equation w.r.to t, we get

`d^2/dt^2theta=d^2/dt^2theta_1=d/dttheta_2` ........................................................... (2)

substituting equations (1) and (2) in equation (A), we get

`d^2/dt^2theta_1+b/m*d/dttheta_1+g/Lsintheta_1`

`d/dttheta_2+b/m*theta_2+g/Lsintheta_1` so

`d/dttheta_2=-b/m*theta_2-g/L*sintheta_1` and `d/dttheta_1 =theta_2`

Converting above two equations in the form of matrix, we get

`d/dt((theta_1),(theta_2))=[[theta_2],[-b/m*theta_2-g/L*sin(theta_1)]]`

In order to run the above function in Python, a list of theta should be made such as

`theta =d/dt((theta_1),(theta_2))=[[theta_2],[-b/m*theta_2-g/L*sin(theta_1)]]`

Objectives:

1. To write a program for second order ODE to show simple harmonic motion.
2. To show damping oscillation of simple pendulum in the form of graphical representation.
3. To animate motion of simple pendulum.

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{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.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 (first Plot represents a single loop output while the second plot represents infinite looped output)

Conclusion:

Second Order differential equations were developed using Python. the necessary animations have also been displayed.