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.And also write a program to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

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 writtendown,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 these systems.

Solution:-

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.

Given,

• L= 1 meter,
• 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.

Coding:

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

#function that returns dtheta/dt 
def function(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 = 1

#initial conditions
theta_0 = [0,3]


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

#solve ode

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

#plot result
plt.plot(t,theta[:,0],'g- ',label=r'$frac{dtheta_l}{dt}=theta2$')
plt.plot(t,theta[:,1],'b--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{l}sintheta_1$')  
plt.xlabel('Plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('Output.png')
plt.show()

#animation
ct = 1

x0 = 0
y0 = 0

for i in theta[:,0]:

        x1=(1*math.sin(i))
        y1=(-1*math.cos(i))

        filename = 'test%05d.png' % ct 
        ct = ct +1

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

Steps used for coding:-

• Import module numpy os np,scipy.integrate import odeint for second ODE,math and for plots matplotlib.
• Define function that returns the value of dtheta/dt.
• Initialize the input values of b,g,m,l.
• Define initial condition form theta i.e., 0 and 3 rad/sec.
• Use np.linspace for creating numeric sequence;np.linspace(start point last point number of values take in betweenstart and end point).
• Solve ODE by use of odeint command.
• Plot result by using pit.plot command to plot graph.
• For animation:- use for loop
• Name the filename
• Use pit.plot command to create every single point and after that stiches allthe pie and create a git

Results: -

Output:- Angular displacement,angular velocity with respect to time.

Pendulam video:

https://drive.google.com/file/d/1KtPXBQ2d9qK7VafDpxT_WxfLDSVpFvBR/view?usp=sharing