Week 3 - Solving second order ODEs

AIM :-

To solve second order differential equation for pendulum by using python program and make animation of 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.

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.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

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

Python program to solve differential equation: -

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

# fuction that return 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    #damping coefficient 
g = 9.81   #gravity in m/s^2
l = 1      #length in meter
m = 1      #mass of ball in kg

# Initial condition 
theta_0 =[0,3]

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

# 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('displacementof ball')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

Algorithms: -

• Import the important libraries required for solving ODE and plotting graphs like matplotlib,odient from scipy.integrate to solve ODE
• Define a function which returns the array dtheta

• Give inputs
• Give initial condition from i.e. 0 to 3rad/s
• Use np.linspace command to create numerical sequence
• Solve ode by using odeint command
• By using plt command plot graphs and give label title to graphs

Plot:-

Python program to simulate pendulum motion:-

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

#function
def model (theta,t,b,g,l,m):
    theta1= theta[0]
    theta2 = theta[1]
    dtheta1_dt = theta2
    dtheta2_dt = -(b/m)*theta2 -(g/1)*math.sin(theta1)
    dtheta_dt = [dtheta1_dt, dtheta2_dt]
    return dtheta_dt

#1nputs

b = 0.05 # damping co-efficient
l = 1 # length of pendulum in m
g = 9.81 # gravity in m/s^2
m = 1 # mass of ball in kg

# initial conditions

theta_x = [0,3]

# time points

t = np.linspace (0,20,150)
#solve ode

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

#plot
X = theta[:,0]
count= 1
for theta in X:
    x1 = 0
    y1 = 0
    x2 = 1*math.sin(theta)
    y2 = -(1*math.cos(theta))
    file ='pendulum%05d.png'%count
    plt.figure()
    plt.plot([-1.5,1.5],[0,0],linewidth=5)
    plt.plot([x1,x2],[y1,y2],linewidth=4)
    plt.plot(x2,y2, 'o' ,markersize =9)
    plt.xlim( [-1.5,1.5])
    plt.ylim([-1.5,1.5])
    plt.savefig(file)
    count = count+1

Algorithms:

• Import the important libraries required for solving ODE and plotting graphs like matplotlib,odient from scipy.integrate to solve ODE
• Define a function with input arguments which returns dtheta_dt
• Give inputs
• Give initial condition from i.e. 0 to 3rad/s
• Use np.linspace command to create numerical sequence
• Solve ode by using odeint command
• By using for loop save image for each count of pendulum
• By using imagemagic in CYWIN command make gif of those images and make animation

Output:

Animation link:-

https://youtu.be/lzusSXkaAuc

Result:-

Angular displacement and angular velocity and it changes over the time due to the effect of damping.