SOLVING 2ND ORDER ODE USING PYTHON

Objective: To solve a second order ODE and create an animation simulating the motion of the pendulum

Problem Statement:

To 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.

To get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum.

Theory:

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.

Similarly, there are hundreds of ODEs that describe key phenomena, and with the ability to solve these equations any system can be simulated.

Differential equations are solved in Python with the Scipy.integrate package using function odeint.

Syntax for odeint function

 y = odeint(model, y0, t)

where,

1. model: Function name that returns derivative values at requested y and t values as dydt = model(y,t)
2. y0: Initial conditions of the differential states
3. t: Time points at which the solution should be reported. Additional internal points are often calculated to maintain accuracy of the solution but are not reported.

Data Given:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2,

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

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b= damping coefficient.

Formulae:

Program:

Code for the main program

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

# 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)

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

#plot results
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
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.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

# Animation of the 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()
   plt.plot([-1.5,1.5],[0,0]) # Fixed End
   plt.plot([x0,x1],[y0,y1]) # String
   plt.plot(x1,y1,'o',markersize = 20) # Pendulum Bob
   plt.xlim([-1.5,1.5])
   plt.ylim([-1.5,1.5])
   plt.savefig(filename)

Procedure:

1. We define the inputs which are passed onto the function later.
2. Then the initial condition for the angle and time is created in the form of variables as mentioned above.
3. Here we are introduced with the odeint function whose syntax is mentioned in the theory, where the parameters such as model, yθ, t (i.e initial condition) are entered.
4. Then the results are plotted for Displacement vs Velocity by using the code mentioned above.
5. This step is a key one because here we try to create an animation for the pendulum, by using a 'for' loop.
6. The code as mentioned above is written and the plotting is done based on the solution of the ode function.
7. The formatting can be done within such as the size of the bob which can be adjusted.
8. Then at each frame the image is captured and stored in the folder by using the savefig() function.
9. Then the animation and the video are created using the additional software which are available on the internet.

Errors:

This error had occurred because the 'm' variable was not passed into the function so the python shouts an error. So by adding the variable 'm' to the argument the error can be overcome.

Result:

This is the plot that describes the variation of displacement and velocity with time.

This is the gif image for the pendulum animation.

Conclusion: We can conclude that with the increase in the time the angular displacement and the angular velocity decreases because of the damping effect which is clearly evident from the animation shown above.

Animation Video Link:

https://www.youtube.com/watch?v=uwZPQGZoDSY