Week 3 - Solving second order ODEs

Aim:

To write a program that solves the given 2 nd order ODE and plot the graphs of solution -motion

The ODE represents the equation of motion of a simple pendulum with damping.

Constants:

L=1 metre,

m=0.1kg

b=0.05.

g=9.81 m/s2.

Initial values:

at t=0 y=[0,3]

program code:

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

theta=[]

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=0.1

fig,ax=plt.subplots()
eq1 = (r"dθ1/dt=θ2")
       #r"\end{eqnarray*}")
ax.text(20, 4, eq1,, fontsize=18,
        horizontalalignment="right", verticalalignment="top")
eq2=(r"dθ2/dt=-(b/m)θ2-(g/l)sinθ1")

ax.text(21, -3.5, eq2,, fontsize=18,
        horizontalalignment="right", verticalalignment="top")

theta_0=[0,3]
t=np.linspace(0,20,150)
theta=odeint(model,theta_0,t,args=(b,g,l,m))
print(theta)

plt.ylim(-4,4)
plt.xlim(0,20)
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}\theta_2-\frac{g}{l}sin\theta_1$')


plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')

plt.show()

The curves as attached below: