Week 3 - Solving second order ODEs

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 to animate the motion of the pendulum
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.

Deliverables :

Submit a report that has the following content:

• Step wise explanation of your code.
• Errors faced while programming.
• How have you overcome those errors.
• Screenshots showing errors thrown in command window.
• Upload an animation on You tube and provide the You tube link for an animation in your report.

Introduction:

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.

Data Given:
L=1 metre, length of the pendulum in m,
m=1 kg, mass of the ball in kg,
b=0.05, damping coefficient.
g=9.81 m/s2, gravitational acceleration in m/s2,

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

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 maintainaccuracy of the solution but are not reported.

Procedure:

1. math and matplotlib.pyplot were imported for using mathematical functions and plotting.
2. numpy was imported for using linspace command
3. odeint was imported to solve the differential equation.

4. The inputs given are as follows:

• damping coefficient (b) = 0.05
• length of the pendulum(l) = 1m
• mass of the pendulum(m) = 1kg

5. Then the initial condition for the angle and time is created in the form of variables as mentioned above. Here we are introduced with the odeint function whose syntax is mentioned in the syntax, where the parameters such as model, yθ, t (i.e initial condition) are entered.The initial angular displacement(0) and velocity(3 rad/sec) were defined by the theta_0 array.

The time input(0 to 20secs) was given by the variable t and was split in 150 parts using linspace command for plotting the curve.

6. A function was defined for calculating the angular displacement and velocity at different time intervals.

The second-orderdifferential equation defining the motion of pendulum

The second order differential equation

((d^2theta)/(dt^2) + (b/m)*(d theta)/dt +(g/l)sin theta=0)

is split into two single order differential equations as

((d theta_1)/ (dt)=theta_2)

((d theta_2)/(dt) = -(b/m)*theta2 -(g/l)*sintheta1)

and defined within the function. odeint was used to solve this differential equation with given input parameters.

7. Here we are introduced with the odeint function whose syntax is mentioned in the theory, where the parameters such asmodel, yθ, t (i.e initial condition) are entered.

8. The displacement and velocity graph against time was plotted using plot command.

9. This step is a key one because here we try to create an animation for the pendulum, by using a 'for' loop.

10.The code as mentioned above is written and the plotting is done based on the solution of the ode function.

11.The formatting can be done within such as the size of the bob which can be adjusted.

11. For creating an animation of pendulum motion, multiple images of different pendulum positions were captured. The different positions were calculated by writing the coordinates of the pendulum as a function of angular displacement. Toconvert the images to animation 'imageMagik' was used.

Program :

#ODE 
import math 
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 
#inputs
b = 0.05 #damping co-eff 
g = 9.81 #gravity 
l = 1    #length of pendulum 
m = 1    #mass of the pendulum 
ct = 1 
theta_0 = [0,3] #initial condition for both ODE of fisrt order here angular displacement-(0) and velocity-(3) 
t = np.linspace (0,20,150) 
#time points or time array where t = 0 sec and integrate it to time = 10 sec
# in between 0-10 we require 150 values
# Function that returns dz/dt
def graph(theta,t): 
	theta1 = theta [0]   #1- 1st order ODE
	theta2 = theta [1]   #2- 1st orfer ODE
	dtheta1_dt = theta2 
	dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1) 
	dtheta_dt = [dtheta1_dt,dtheta2_dt] 
	return dtheta_dt
theta = odeint (graph,theta_0,t) #solving ODE
plt.plot(t,theta[:,0], 'b-',label= r'$\frac{d\theta1}{dt} =\theta2$') 
plt.plot(t,theta[:,1], 'r--', label= r'$\frac{d\theta2}{dt} =-\frac{b}{m}\theta2\frac{g}{L}sin\theta1$')
plt.legend()
plt.ylabel('Displacement and Velocity') 
plt.xlabel('Time in Sec')
plt.legend(loc='best')
plt.show()
for i in theta[:,0]: #loop for pendulum animation 
	X1 = l*math.sin(i) 
	Y1 = -l*math.cos(i)
	plt.figure() 
	plt.xlim([-1,1]) 
	plt.ylim([-1.5,0.1])
	plt.ylabel('Plot') 
	plt.xlabel('Time in Sec')
	plt.grid()
	plt.plot([0,X1],[0,Y1]) 
	plt.plot([-0.5,0.5],[0,0]) 
	plt.plot(X1,Y1,marker='o',markersize=20) 
	plt.savefig(str(ct)+'.png') 
	ct=ct+1

12. Then the animation and the video are created using the additional software which are available on the internet.

Error :

This error had occurred because the handle with label is not given. It is not neccesary to run the code. It can be overcomed by giving the labels to legend.

Results:

1. The graph below shows the variation of angular displacement and velocity (on y-axis) with respect to time (on x-axis).
2. The displacement curve (blue line) can be seen starting at value 0 as specified in the given parameters. Similarly, thevelocity curve (orange line) starts at 3 as given in the input.
3. The amplitude can be seen decreasing in both the graph and the gif. This occurs due to the damping coefficient given inthe input parameters.

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.