Week 3 - Solving second order ODEs

Objective:

To solve second order differential equation

Description:

In engineering, ODE is used to describe the transient behaviour of a system

Example: Simple pendulum

A simple pendulum is a bob of mass (m) suspending to string of length (l) from a pivot so that it can swing freely, when a bob is displaced sideways from its equilibrium position, it is subjected to a restoring force due to its acceleration due to gravity(g) ,it will accelerate backward toward the equilibrium position. when released the restoring force acting on mass causes to oscillate about its equilibrium position, swinging back and forth. The time for one complete swing, a right swing & left swing is called the period. The time period depends on the length of the pendulum and also to as slight degree on amplitude.

Systematic representation of a simple pendulum

Time period (T)=2*π*`sqrt(l/g)`

l= length of the string

g=acceleration dude to gravity

Second order differential equation is given by:

`theta`= angle made by th string

m= mass of the bob

b= damping coefficient

g= accelertaion due to gravity

l=length of the string

Code explanastion:

Inputs:

L=1 meter,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Time =0 to 20sec

Angular displacement=0

Angular velocity =3 rad/sec

1.Import math,matplotlib,numpy,scipy,imageio modules

2.Define a second order differential function

3.Convert second order ODE in first order ODE

4.import odeint from scify module

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

Syntax: y=odient(model,y0,t,args)

model:: Function name that returns derivative values at requested y and t values as

dydt = model(y,t)

y0: Initial conditions of the differential states

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.

args: An optional fourth input is args that allows additional information to be passed into the model function. The args input is a tuple sequence of values. The argument k is now an input to the model function by including an addition argument

6.linspace command is in numpy module is used to extract the time points between start time and end time

7.Plot the graph Angular displacement vs time and Angular velocity vs time

8.xlabel and ylabel is used to rename the axis

9.axis() is used to scale the axis

10.append() command is used to add the list of values into an empty set

11.imread() command is used to read the list of images

12.imsave() command is used to save the list of images into a required format. In this code .gif format used

Code for plotting Angular dispalcement Vs Time and Angular velocity vs time

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

# function that retuns 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

b=0.05
g=9.81
l=1
m=0.1

#initial conditions
theta_0=[0,3]	
# time points

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

#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\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin|theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

Code for plotting Animation of a simple pendulum as a function of time:

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

# function that retuns 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  #Damping coefficient
g=9.81   #accelertaion due to gravity in m/s^2
l=1    #length of the string in metres
m=0.1   #mass of the bob in kg

#initial conditions
theta_0=[0,3]	
# time points

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

#solve ODE

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


# Animation of a simple pendulum
ct=1 #contor variable
x0=0
y0=0

x=theta[:,0]
images=[]

for i in x:
	x1=x0+l*math.sin(i)
	y1=-(y0+l*math.cos(i))

	filename=str(ct)+'.png'

	plt.figure()
	plt.plot([-1.5,1.5],[0,0])
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='o',markersize='12')
	plt.axis([-1.5,1.5,-1.5,0.5])
	plt.savefig(filename)

	images.append(s.imread(filename)) 
	ct=ct+1
	s.mimsave('simple_pendulum.gif',images)

Results:

1.Angular displacement vs time and Angular velocity vs time

2.Animation of a Simple pendulum as a function of time