Week 3 - Solving second order ODEs

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

Problem Statement:

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.

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.

• For each value in time
• Create an image that represents the position of the pendulum

Solution Procedure:

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.

First of all convert second order ode into the two second order ode.

Then assume theta as a solution vector

So theta1=angular displacement

And theta2=angular velocity

Import libraries:

Especially scipy module

function that gives or returns dtheta/dt

Defining the values for damping co.eff ,gravity,length of pendulum,mass of ball

initial conditions angular displacement =0 and angular velocity=3

time points 0 to 20 sec bu using numpy module and using linspace command

solve ode by odeint

plot result

for animation take counter variable and save file name by for loops

ploting and getting images for different angle values

python code is following:

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

#function that gives or returns dtheta/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    

b = 0.05 #damping co.eff 
g = 9.81  #gravity
l = 1  #length of pendulum
m = 1  #mass of ball

#initial conditions  angular displacement =0 and angular velocity=3
theta_0 = [0,3]   

#time points  0 to 20 sec
t = np.linspace(0,20,150)    

#solve ode    
theta = odeint(model,theta_0,t,args = (b,g,l,m))    

#plot result  
plt.plot(t,theta[:,0],'g-',label=r'Displacement')  
plt.plot(t,theta[:,1],'r--',label=r'Velocity')  
plt.xlabel('Plot')  
plt.ylabel('Time')  
plt.legend()  
plt.show()    

#animation  
ct = 1    

x0=0  
y0=0    

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

#ploting and getting images for different angle values 

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

Results and Discussion:

Error and output:

Animation:

https://youtu.be/YrwLJIgtC6M