Challenge

Week 3 - Solving second order ODEs

https://youtu.be/8G1IV2tbSCc

"""
Objective - To solve the Second order ordinary Differential Equation. the taken equation is for simple pendulum.
Hence derive trace the path of simple pensulum
Create Animation for Pendulum
"""
#importing the operator from Library
import numpy as np
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt 
import matplotlib.animation as animation




# define the function and then  returning its value
# this fuction is second order differential equation for simple pendulum.
# this is asper one given for exercise
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.02  # damping coefficient 
g = 9.81 # gravity (m/s2)
l = 1 #length of string (m)
m = 0.1 #mass of ball (Kg)

# defining time; this is the total instant thru which the theta value is observed
t = np.linspace(0,20,100)

# defining initial angle condition
# where math.pi/2 is position and 0 is the velocity, thus defining the initial condition
theta_0 = [math.pi/2,0]

# defining the ode
'''
this odeint is the inbuilt fuction in the python for solving differential equations
the 'model' is defined function for second order derivative, theta_0 initial condition, t the time instances limiting the upper and lower limit for derivation
args - are the paramters that will be substituted from the consideration.

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


# plot results
# below plots give us the sinosidal waves with respecte to TIME,
# theta[:,0] =  velocity  - displacement WRT Time i.e first derivative
#theta[:,1] = acceleration - derivative of velocity with repsect to TIME i.e seconf derivative of displacement wrt time
plt.plot(t, theta[:,0], 'b-', label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t, theta[:,1], 'r--', label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('time')
plt.legend(loc='best')  # defining the position for the legends
plt.show()


# Now below are the CODE for creating Animation
# animation for the pendulum

image = plt.figure(figsize=(5,5), facecolor = 'w') # this creates main image
ar = image.add_subplot(1,1,1) #this creates the subplot with the size 1x1x1
wall, = ar.plot([-1,1],[0,0], color = 'black', lw = 5) # wallof pendulum
plt.rcParams['font.size'] = 20


#creating an array for appending the calculated value from forloop into it
string = []



"""
as theta value is obtainedformthe function odeint, for 300instance. This is becasue we have considered 300 steos of TIME t
hence using these theta instant value to plot pendulum wrt to x&y axis cordinate
"""
for i in range(len(theta)):

	#plotting the string with x and y co-ordinate generated. as we are subplotting it,we need to make sure that ',' is used after name as syntax requirement
	string_tmp, = ar.plot([0,np.sin(theta[i,0])], [0, -np.cos(theta[i,0])], color = 'k', lw = 2)

	#plotting markker - plots end point ofthe string
	end_mark, = ar.plot(np.sin(theta[i,0]), -np.cos(theta[i,0]), 'o',markersize = 10, color = 'r')

	#plotting time simulation
	time = ar.text(-0.9,0.25, 'Time = %.1fs' % t[i])

	# As we have taken special empty array,in that we are going put all the appended value i order to store it
	string.append([string_tmp, end_mark, time, wall])
	
# this syntax sets the limit  for thex data and y data
ar.set_aspect('equal', 'datalim')


"""
there are numeraous way of creating annimation like capturing the images and then stiching it.
this will require another tool to stich and create the .gif
IN this excercise we used different method. 
for 'ffmpeg'is pip installed in the system as a writer
also we need to import animation from the MATPLOTLIB libraray i.e 'artistanimation'
"""
#saving animation
video = animation.ArtistAnimation(image,string,interval = 50) # this plots image and data of 'string' that is captured
fn = 'wall_Pendulum_animation'
video.save(fn + '.mp4', writer = 'ffmpeg', fps=1000/50) # this is general syntax that must be followed


Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.