2022-01-31-135759.py

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

M=5
m=2
g=9.81

def Atwood(y,t,M,m,g):

    #theta'(t)=omega(t)
    #radius'(t)=Rho(t)

    #omega'(t)=(1/radius(t))(-g*sin(theta(t))-2(Rho(t)*omega(t))
    #Rho'(t)=(1/(M+m))(m*radius(t)*(omega(t)^2)-M*g+m*g*cos(theta(t)))

    theta, omega, radius, Rho = y

    dydt = [omega, (1/radius)*(-g*(np.sin(theta))-2*(Rho*omega)), Rho, (1/(M+m))*(m*radius*(omega**2)-M*g+m*g*(np.cos(theta)))]

    return dydt

#y0=[theta, change in theta, radius, change in radius]
y0 = [0 ,2 ,1 ,1]
#time points
t = np.linspace(0,10,101)
#Solving the ODE
sul = odeint(Atwood, y0, t, args=(M,m,g))

plt.plot(sul[:,0], sul[:,1], 'r', label='theta (x) vs theta dot (y)') #Phase plot
plt.legend(loc='best')
plt.grid()
plt.show()