code

import math
import matplotlib.pyplot as plt

def Df(F_D,LIM):
    w=[0 for i in range(0,LIM)]
    theta=[0 for i in range(0,LIM)]
    t=[0 for i in range(0,LIM)]
    w[0]=0;theta[0]=0.2;Omega_D = 2/3;dt=0.04;g=9.8;l=9.8;q=0.5

    for i in range(0,LIM-1):
        w[i+1]=w[i] + (-(g/l) * math.sin(theta[i]) - q*w[i] + F_D * math.sin(Omega_D * t[i]) )*dt
        theta[i+1]=theta[i]+w[i+1]*dt
        if theta[i+1] > math.pi:
            theta[i+1]=theta[i+1]-2*math.pi
        elif theta[i+1] < -math.pi:
            theta[i+1]=theta[i+1]+2*math.pi
        else:
            theta[i+1]=theta[i+1]
        t[i+1]=t[i]+dt
    plt.subplot(221)    
    plt.plot(t,w)
    plt.subplot(222)
    plt.plot(t,theta)
    plt.subplot(212)
    plt.scatter(theta,w,s=6)
    plt.show()
Df(1.2,1500)