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 :-

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When we do this we will get the position of the pendulum as a function of time. We can use this information to animate the motion of the pendulum, as shown in the below figure.

Code :-

# To solve 2nd order ODE using python

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

#function that returns dtheta/dt

def function(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.5
g=9.81
l=1
m=1

# initial conditions
theta_0 = [0,3]

# input for time
t = np.linspace(0,20,200)

theta = odeint(function,theta_0,t,args=(b,g,l,m))

#plot result
plt.plot(t,theta[:,0],'-g',label = r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'b--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{1}sin\theta_1$')
plt.xlabel('Plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('output.png')
plt.show()

# To create Animation
ct = 1
x0 = 0 
y0 = 0

for i in theta[:,0]:
	x1=x0+(1*math.sin(i))
	y1=y0+(-1*math.cos(i))

	filename = str(ct)+'.png'
	ct = ct+1

	plt.figure()
	plt.plot([-1,1],[0,0],"black")
	plt.plot([x0,x1],[y0,y1],"red")
	plt.plot(x1,y1,'o',markersize = '20')
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.0])
	plt.savefig(filename)

Solution Procedure :-

First of all, we import the libraries that are going to be used in this program. ‘math’ module for basic math commands, ‘matplotlib.pyplot’ module for plotting purposes, ‘numpy’ module for arrays and ‘odeint’ module from ‘scipy.integrate’ for integrating our second order ODE.

Before assigning the input values we write a function named ‘model’ with inputs as ‘theta’, ‘time’, ‘b’, ‘g’, ‘l’ and ‘m’. This function breaks our second order ODE into 2 first order ODEs and returns the value of dtheta_dt.

Then we assign the input values of damping coefficient (b), gravitational acceleration (g), length of pendulum (L) and mass of the bob (m).

After that we define a variable theta_0 which contains the initial angular displacement = pi/2 and initial angular velocity = 0. After that we assign the value of time span as 0 to 20 sec by the help of linspace command with 200 intervals.

Now to solve the ODE, we use the command ‘odeint’ which is going to integrate the ODEs and inside the command, we call the function ‘model’ with all the inputs and the arguments.

Then we use the plot command to generate the plot between displacement and time and velocity and time.

Now to create an animation of pendulum, first we define the initial values of the coordinates, which is going to be a fixed point. Then we used a ‘for’ loop with loop counter variable (i) with size equal to the length of time (t) array. Also count (ct) variable is also addressed to create separate file for each case.

Then we define the position of the end point of the string. To create the mass at one end of the pendulum, we use ‘o’ inside the plot command for end point of string, then we choose the color and choose the size by using ‘markersize’ command.

To name the file in a definite order, we create a variable named ‘filename’ with str of ct and extension as ‘.png’ and use ‘savefig’ command to save the files.

Now when we run the code, all 200 images are saved in the working folder. We stitched these images together to create the required GIF file.

Result :-

1. Displacement and velocity plot with respect to time :

2. GIF of animation of Simple Pendulum :

3. Onedrive link for the animation video :

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

Errors -

1.Error 1

The Error is solved by giving the value "dtheta2_dt = -(b/m)*theta2-(g/l)*math.sin(theta1)" in the line 14.

2.Error2

The Error is soved by giving the value "x1=x0+(1*math.sin(i))" in the line 47.

Conclusion :-

Finally, the second order ODE is solved using Python and the animation for a simple pendulum is created. By simply changing the inputs used, we can understand the effect of different parameters.