Solving second order ODEs using python

AIM:

To solve 2nd order ODE's.

OBJECTIVES:

1. To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

2. To animate the motion of the pendulum.

PROBLEM STATEMENT:

Write a program in Python that will simulate the pendulum motion.

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 you do this you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum.

THEORY:

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.

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.

PROGRAM:

"""
Program for solving 2nd order ODE
Program by Saurabh Kshirsagar
"""
import numpy as np 
import math
import matplotlib.pyplot as plt 
from scipy.integrate import odeint

#Function that returns dz/dt
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 condition
theta_0 = [0,3]

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

#Solving ODE
theta = odeint(model,theta_0,t,args=(b,g,l,m))

#Plotting results
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.xlabel('Time')
plt.ylabel('Plot')
plt.legend(loc='best')
plt.show()

#Animation of Pendulum
ct = 1
x0 = 0
y0 = 0

for i in theta[:,0]:
	x1 = l*math.sin(i)
	y1 = -l*math.cos(i)
	filename = str(ct) + '.png'
	ct = ct+1

	plt.figure()
	#fixed end
	plt.plot([-1.5,1.5],[0,0])
	#string
	plt.plot([x0,x1],[y0,y1])
	#bob
	plt.plot(x1,y1,'o',markersize=22)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)

EXPLANATION:

1. In this program we imported numpy to work with arrays, math to solve some trignometric functions, matplotlib to plot our data and odeint from scipy to solve ODEs.
2. We defined a function under which 2nd order ode is converted to 1st order ode.
3. Inputs, initial conditions and time points were given as per the problem definition.
4. We used odeint which was imported from scipy to solve our ODE.
5. Results are plotted for angular displacement and angular velocity vs time using matplotlib module.
6. Finally we created the animation using magicks software.

ERRORS FACED:

OUTPUTS:-

1. Plot of variation of agular velocity and angular displacement with time.

2. Animation of motion of pendulum.

CONCLUSION:

From the plot and animation we can see that with increase in time the angular velocity and angular displacement decreases.

Animation Youtube Link:- https://youtu.be/lCiGQATZKj0