from numpy import * from matplotlib.pyplot import * # TASK 1.2 class CoupledPendulum: """A class that simulates a spring pendulum.""" def __init__(self, dt, k1, k2, kc, x1_init, x2_init, v1_init = 0.0, v2_init = 0.0): # spring constant self.k1 = k1 self.k2 = k2 self.kc = kc # time step self.dt = dt # time self.t = 0.0 # current position self.x1 = x1_init self.x2 = x2_init # current velocity self.v1 = v1_init self.v2 = v2_init def step(self): """Perform a single time step.""" # compute forces F1 = -self.k1*self.x1 - self.kc*(self.x1-self.x2) F2 = -self.k2*self.x2 - self.kc*(self.x2-self.x1) # velocity step self.v1 += F1*self.dt self.v2 += F2*self.dt # compute new position self.x1 += self.v1*self.dt self.x2 += self.v2*self.dt # increase time self.t += self.dt def getEnergy(self): """Compute and return the energy components.""" # kinetic energy e1 = 0.5*self.v1**2 e2 = 0.5*self.v2**2 # potential energy e1 += 0.5*self.k1*self.x1**2 e2 += 0.5*self.k2*self.x2**2 ec = 0.5*self.kc*(self.x1-self.x2)**2 # total energy e = e1+e2+ec return e, e1, e2 def simulate(self, maxTime): """Simulate the pendulum up to maxTime. Returns the measurements as a NumPy array.""" measurements = [] while self.t < maxTime: self.step() e, e1, e2 = self.getEnergy() measurements.append((self.t, self.x1, self.x2, e, e1, e2)) return array(measurements) p = CoupledPendulum(0.1, 1.0, 0.7, 0.2, 0.3, 0.2) m = p.simulate(30.0) # Create an empty figure figure() # Plot of the positions vs. time subplot(311) xlabel("Time t") ylabel("Position x") plot(m[:,0], m[:,1], 'o-') plot(m[:,0], m[:,2], 'o-') # Plot of the energies vs. time subplot(312) xlabel("Time T") ylabel("Energy E") plot(m[:,0], m[:,4], 'o-') plot(m[:,0], m[:,5], 'o-') subplot(313) xlabel("Time T") ylabel("Total Energy E") plot(m[:,0], m[:,3], 'o-') # Task 1.2.2 class CoupledPendulum2 (CoupledPendulum): def step(self): """Perform a single time step.""" # compute forces F1 = -self.k1*self.x1 - self.kc*(self.x1-self.x2) F2 = -self.k2*self.x2 - self.kc*(self.x2-self.x1) # compute new position self.x1 += self.v1*self.dt self.x2 += self.v2*self.dt # velocity step self.v1 += F1*self.dt self.v2 += F2*self.dt # increase time self.t += self.dt p2 = CoupledPendulum2(0.1, 1.0, 0.7, 0.2, 0.3, 0.2) m2 = p2.simulate(30.0) figure() xlabel("Time T") ylabel("Total Energy E") plot(m[:,0], m[:,3], 'o-') plot(m2[:,0], m2[:,3], 'o-') # Now show the graphs show()