Week 3 - Solving second order ODEs

TITLE: TO SOLVE 2ND ORDER DIFFERENTIAL EQUATION USING PYTHON

OBJECTIVE: To solve 2nd ODE and create animation of pendulum.

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.

Formula used:

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.

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Steps of coading:

• At first i import the libraries.
• Then I define the functions of theta,L,b,m and g.
• Then I use the above formula to the get the graph.
• Then I use plt command to plot the graphs.
• Then I write code to get the plots for pendulum.
• Here I use for loop to get 150 diffrent graphs.
• Then I use 'test%05d.png' % ct' command to assemble all the graphs in order.
• Then I use image magick to stitch the graphs and get an animation file.

CODE:

#To slolve 2nd order diffrential equation import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt import math #function that returns dz/dt def model(theta,t,b,g,l,m): theta1 = theta[0] theta2 = theta[1] dtheta1_dt = theta2 dteheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1) dtheta_dt = [dtheta1_dt ,dteheta2_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,10,150) #solve ODE theta = odeint(model,theta_0,t,args = (b,g,l,m)) print(theta) #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}\theta2 -\frac{g}{l}sin\theta_1$') plt.ylabel('plot') plt.xlabel('time') plt.legend() plt.show() #program for the animation of pendullum theta1= theta[:,0] ct = 1 x1 = 0 y1 = 0 for t1 in theta1: #Angular velocity x2 = l*math.sin(t1) y2 = -l*math.cos(t1) filename = str(ct) + '.png' filename = 'test%05d.png' % ct ct = ct + 1 plt.figure() plt.plot([-2,2],[0,0]) plt.plot([x1,x2],[y1,y2]) plt.plot(x2,y2,'o',markersize = 10) plt.xlim([-2,2]) plt.ylim([-2,1]) plt.savefig(filename)

RESULTS: