from __future__ import division from math import * from pylab import * # define a function that given amplitude of force, gives angle,avelo and t def change_amp(F): # define some constants q = 0.5 l = 9.8 g = 9.8 f = float(2/3) dt = 0.04 theta0 = 0.2 omega0 = 0 angle = [theta0] avelo = [omega0] t = [0] # move the pendulum while t[-1] < 60: avelo_new = avelo[-1] - ((g/l)*sin(angle[-1]) + q*avelo[-1] - F*sin(f*t[-1]))*dt avelo.append(avelo_new) angle_new = angle[-1] + avelo[-1]*dt while angle_new > pi: angle_new -= 2*pi while angle_new < -pi: angle_new += 2*pi angle.append(angle_new) t_new = t[-1] + dt t.append(t_new) return angle,avelo,t angle_0 = change_amp(0)[0] avelo_0 = change_amp(0)[1] t_0 = change_amp(0)[2] angle_1 = change_amp(0.5)[0] avelo_1 = change_amp(0.5)[1] t_1 = change_amp(0.5)[2] angle_2 = change_amp(1.2)[0] avelo_2 = change_amp(1.2)[1] t_2 = change_amp(1.2)[2] plot(t_0,angle_0,label='FD=0') plot(t_1,angle_1,label='FD=0.5') plot(t_2,angle_2,label='FD=1.2') legend(loc='upper right') title('angle versus time') xlabel('time(s)') ylabel('theta(radians)') xlim(0,60) grid(True) show()