simplependulumphaseplot.py

#Import recquired models
import numpy as np
import matplotlib.pyplot as plt

#This allows us to code for ODEs
from scipy.integrate import odeint

#Setting values
g = 9.81
L = 2

#Function- We convert the SHM ODE of a pendulum into 2 first order ODEs
def model(y,t,g,L):
    
    #theta'(t) = omega(t)
    #omega'(t) = -(g/L)*(sin(theta(t)))
    
    theta, omega = y
    dydt = [omega, -(g/L)*(np.sin(theta))]
    return dydt
    

#initial conditions
#sets that at time=0, angular velocity is 3 radians per second
y0 = [0-0.1,0.0]

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

#solve ODE
sol = odeint(model, y0, t, args=(g,L))

#plotting results
plt.plot(sol[:,0], sol[:, 1], 'b')
plt.ylabel('Angular velocity (radians per seconds)')
plt.xlabel('Angular displacement (radians)')
plt.title('Fig 3')
plt.grid()
plt.show()
plt.show()