Week 3 - Solving second order ODEs

#second order ODE
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#function that returns 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.02
g=9.81
l=0.1
m=0.1

#initial condition
theta_0=[0,5]

#time points
t=np.linspace(0,10,150)

#solve ODE
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#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.xlabel('time')
plt.ylabel('plot')
plt.legend(loc='best')
plt.show()

#counter variable 
ct=1
for theta in theta[:,0]

#coordinates of fixed points of string
x0=0
y0=0
#coordinats of fixed points of string
x1=1*math.sin(theta)
x2=-1*math.cos(theta)

filename=str(ct)+'png'
ct=ct+1

plt.figure()
#plotting horizontal line
plt.plot([-0.5,0.5],[0,0],'g',linewidth=5)
#plotting string of pendulum
plt.plot([x0,x1],[y0,y1],'r')
#plotting ball of pendulum
plt.plot(x1,y1,'o',markersize=30)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.title('animation of simple pendulum using python')
plt.savefig(filename)

OUTPUT:-

link for simulation video of simple pendulum:-

https://youtu.be/HU3KfSueBT4

explanation of code:-

1. Now we are solving the ODE as mentioned above and we split that ode into two single order ODE’s one for velocity and one for displacement.
2. So, when we need to solve the first ODE, we need to provide the initial conditions for the displacement and velocities
3. Initial conditions which are displacement and velocity are given in the variable 'theta_0'. These conditions are further used to solve the ODE
4. The next step includes giving the time points for the ODE. These points are used for integrating while solving the ODE. The 'linspace' command is used to divide the time interval into required parts
5. Plotting the 'Velocity & displacement vs Time' graph. The 'plot' command is used to plot the graph.
6. To create an animation for the pendulum a conter variable 'ct' which is equaled to 1 is created.
7. math module is imported in order to perform trignometric functions in the program
8. numpy module is imported in order to perform scientific programming in the program
9. odeint module is imported in order to perform ordinary differential equations in the program
10. A function is defined as model which is function of t and all input variable
11. savefig() command is used to save the frames with already defined filename
12. plt.figure() opens new frame window for each frame
13. plt.plot() is used to plot the position of string
14. plt.xlim and plt.ylim is used to define the limiting values of both axis of graphs