Solving second order ODEs
Objective: Determine the problem and solve it though code
Problem statement: Defined in the following steps with solutions.
Step wise explainiation has given before every code. And python file is been attached, for reference.
By useing FCP- Final Cut Pro, I have stitched the images and created video of containin 13 seconds, which is uploaded in youtuble. And link has been provinded in the end of the document.
While scripting the code, came across many errors and solved it. I have attached screenshot of the errors.
Errors are like, human errors.
INSERT SAMPLE CODE
#Program to solve second order ODEs and simulate the motion of the pendulum ball
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#Function that return dz/dt
def model(theta, t, b, g, l, m):
theta_1 = theta[0]
theta_2 = theta[1]
d_theta1_dt = theta_2
d_theta2_dt = -(b/m)*theta_2 - (g/l)*math.sin(theta_1)
d_theta_dt = [d_theta1_dt, d_theta2_dt]
return d_theta_dt
b = 0.5 #Damping coefficient
l = 1 #Length of the pendulum in metres (m)
g = 9.81 #Gravity in (m/s^2)
m = 1 #Mass of the ball in kilograms (Kg)
#Initial conditions
theta_0 = [0 ,5]
#Time points
t = np.linspace(0,10,150)
#Solve ODEs
theta = odeint(model, theta_0,t,args = (b, g, l, m))
# Plotting the graph
plt.plot(t, theta[:,0],'b-',label=r'$frac{dtheta_1}{dt}=theta_2$')
plt.plot(t,theta[:,1],'r--',label=r'$frac{dtheta_2}{dt}=-frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.xlabel('Time (second)')
plt.ylabel('Displacement')
plt.title('Harmonic motion of Simple Pendulum')
plt.grid(axis='y')
plt.legend(loc = 'best')
plt.show()
#Now we are going to code for animating ball and string. Along with frames of motion of ball movement(damping)
ct = 1
for i in theta[:,0]:
x1 = 0
y1 = 0
x2 = x1+(l*math.sin(i))
y2 = y1+(-l*math.cos(i))
Damping = str(ct)+'.png'
ct = ct+1
plt.figure()
plt.plot([x1,x2],[y1,y2],linewidth = 2) #string of the pendulum
plt.plot(x2,y2, 'o', markersize = 20) #Bob of the pendulum
plt.plot([-1,1],[0,0], linewidth = 5) #Rigid support
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.0])
plt.title('Damping Pendulum motion')
plt.savefig(Damping)
Animation video of damping pendulum
Youtube video link: https://youtu.be/MxDL9LWo2tI