Week 3 - Solving second order ODEs

2^nd order ODE by using Python:

Aim :

To write a program in Python that will simulate the pendulum motion.

Given the governing equation :

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.

Given data :

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Objective:

• To write a program that solves the given ODE.
• To write a program in Python that will simulate the pendulum motion.
• To animate the motion of the pendulum by using the given information.

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 Newton’s 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.

To solve the given 2^nd order ODE, we, first, must break down that equation in a single order ODE.

We will use these two ODE to calculate the required values.

Before writing the program, we must import Scipy library to calculate the integration part. Otherwise, your program is not going to execute.

First, we will write down the given input information such as g (gravity in m/s2), L( length of the pendulum in m), m( mass of the ball in kg), b(damping coefficient).

Afterward, we will write down the given initial condition such as angular displacement=0, angular velocity=3 rad/sec. and time is given for the motion such as 0-20 sec.

We will create a function as a model to execute and return the value of the given two derivative equation to the odeint including initial conditions and time points where it performs the calculation and an optional fourth input is args that allows additional information to be passed into the model function.

After the odeint, we will plot the values of angular displacement and angular velocity which will odeint provide against the time points respectively.

Lastly, we will create the frames to perform animation of the pendulum by using the given information i.e. angular displacement where we used the ‘for ‘loop.

The following are the code of 2^nd order ODE, graph, and animation of a pendulum by using Python:

Code:

"""
2nd Order ODE using python

BY Ajinkya Nikam
"""

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

#Function which will return the value of dz/dt

def model(Theta, t, b, g, l, m):
	 Theta1 = Theta[0]       
	 Theta2 = Theta[1]
	 dtheta1_dt = Theta2                                  # First ODE
	 dtheta2_dt = -(b/m)*Theta2 - (g/l)*math.sin(Theta1)  # Second ODE
	 dtheta_dt =[dtheta1_dt, dtheta2_dt]
	 return dtheta_dt


# Given values
b =0.05   # damping coeficient
g =9.81   # gravity in m/s2,
l =1      # length of the pendulum in m
m =1      # mass of the ball in kg


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

#initial condition
Theta_rest = [0, 3] 

# Time (evenly spaced numbers over a specified interval)

t = np.linspace(0, 20, 150)

# Solution of ODE

Theta = odeint(model, Theta_rest, t, args = (b, g, l, m))   

#Plotting of graph  

plt.plot(t, Theta[:,0],'b', label = r'$\frac{d\theta_1}{dt}=\theta2$')   # It gives the angular displacement values from the first ODE
plt.plot(t, Theta[:,1],'r--', label = r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$') # It gives the angular velocity values from the second ODE
plt.legend(loc = 'best')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.show()

# To create frames for animation

Ct = 1
for i in Theta[:,0]:
	X0 = 0
	Y0 = 0
	X1 = -l*math.sin(i)
	Y1 = -l*math.cos(i)
	Ct = Ct + 1
	Filename = str(Ct) + '.png'

	plt.figure() 
	plt.plot([-3, 3],[0,0], linewidth = 10, color = "k")
	plt.plot([X0, X1],[Y0, Y1], linewidth = 2, color = "y")
	plt.plot([X0, X1],[Y0, Y1],'D', color = 'w')
	plt.plot(X1,Y1,'o', markersize = 25 , color = 'm')
	plt.savefig(Filename)
	plt.xlim([-2, 2])
	plt.ylim([-2, 2])

Graph:

Animation:

Click here

Conclusion:

We write down the program to solve the 2^nd order ODE and to simulate a pendulum by using python.