from numpy import * from matplotlib.pyplot import * # TASK 1.1.1 class Pendulum: """A class that simulates a spring pendulum.""" def __init__(self, k, dt, x_init = 0.0, v_init = 0.0): # spring constant self.k = k # time step self.dt = dt # time self.t = 0.0 # current position self.x = x_init # current velocity self.v = v_init def step(self): """Perform a single time step.""" # compute force F = -self.k*self.x # compute new velocity self.v += F*self.dt # compute new position self.x += self.v*self.dt # increase time self.t += self.dt def getEnergy(self): """Compute and return the energy components.""" # kinetic energy ek = 0.5*self.v*self.v # potential energy ep = 0.5*self.k*self.x**2 # total energy e = ek+ep return e, ek, ep 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, ek, ep = self.getEnergy() measurements.append((self.t, self.x, e, ek, ep)) return array(measurements) p1 = Pendulum(1.0, 0.1, 0.3) m1 = p1.simulate(10.0) p2 = Pendulum(0.7, 0.1, 0.1) m2 = p2.simulate(10.0) # Create an empty figure figure() # TASK 1.1.2 # Plot of the positions vs. time subplot(221) xlabel("Time t") ylabel("Position x") plot(m1[:,0], m1[:,1], 'o-') plot(m2[:,0], m2[:,1], 'o-') # Plot of the energies vs. time subplot(222) xlabel("Time T") ylabel("Total Energy E") plot(m1[:,0], m1[:,2], 'o-') plot(m2[:,0], m2[:,2], 'o-') subplot(223) xlabel("Time T") ylabel("Kinetic Energy E") plot(m1[:,0], m1[:,3], 'o-') plot(m2[:,0], m2[:,3], 'o-') subplot(224) xlabel("Time T") ylabel("Potential Energy E") plot(m1[:,0], m1[:,4], 'o-') plot(m2[:,0], m2[:,4], 'o-') # TASK 1.1.3 class PendulumForwardEuler(Pendulum): def step(self): """Perform a single time step.""" # compute force F = -self.k*self.x # compute new position self.x += self.v*self.dt # compute new velocity self.v += F*self.dt # increase time self.t += self.dt figure() p1 = Pendulum(1.0, 0.1, 0.3) m1 = p1.simulate(30.0) p2 = PendulumForwardEuler(1.0, 0.1, 0.3) m2 = p2.simulate(30.0) plot(m1[:,0], m1[:,2]) plot(m2[:,0], m2[:,2]) # Now show the graphs show()