Week 3 - Solving second order ODEs

Objective:

The main objective of this challenge is to write a python code for the 2nd order differential equation and plot the graph between time and ang velocity and also create a animation of pendulum.

Problem Statement:

 Given 2nd order differential equation:

`[d^2theta]/dt^2 +b/m [d theta]/dt +g/l*sin theta = 0`

Converting it to multiple first order equations,

`theta = theta1`------eq1

`[d theta]/dt =[d theta1]/dt = theta2`------eq2

`[d^2theta]/dt^2 = [d^2theta1]/dt^2 =d/dt([d theta1]/dt) = d/dt(theta 2)` -------eq3

from all the above equations we simplified the equation as below;

`[d theta2]/dt +b/m theta 2 +g/l*sin theta 1 = 0`

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.

t=0 -20sec

using the above equation the problem need to be solved.

Solution Procedure:

1.  import required libraries including scipy module
2. create a function using the input parameters so that it simplify the code
3. give the input parameters
4. use linspace command to divide the values from 0 to 20 in 300 types
5. now solve the ode
6. plot the results on graph and name them on x and y axis
7. now use the basic math equation to find out the theta values of the pendulum
8. create loop so that it calculates the values of theta in all the cases
9. plot the results on graph
10. use save command to save all the 300 files of different equations
11. use imagemagic to attach them into one to create a smooth gif.

Results and Discussion:

 using the above equations and values we get the required graph between time and angular velocity and also pendulum animation.

Errors faced:

1. indentation error 
2. type error
3. run time error

 

 

Python code:

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

#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=1
ct=1
# initial condition
theta_0 = [0,5]
#time points 
t= np.linspace(0,20,300)
#solve ode
theta =odeint(model,theta_0,t,args = (b,g,l,m))
#plot results
plt.plot(t,theta[:,1],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,0],'r--', label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.legend(loc='best')
plt.xlabel('time')
plt.ylabel('angular velocity')
plt.show()

a = theta[:,0]
b = theta[:,1]
ct = 1

for i in a:
	x0 = 0
	y0= 0
	x1 = 1*(math.sin(i))
	y1 = -1*(math.cos(i))
	filename = str(ct) + '.png'
	ct=ct+1
	plt.figure()
	plt.plot([-0.5,0.5],[0,0])
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='.',markersize=35)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)

Animation of the pendulum:

https://drive.google.com/file/d/1KJfigXx_W0mlM47ts0SM3ev8fdrNK-y3/view?usp=sharing