1.  The mathematical pendulum (20 points)

In class, we have developed a code to simulate the motion of a pendulum. Use this code to investigate the following questions:

(a)  (5 pts) We want to quantify the error of the Euler and the velocity Verlet algorithm in more

detail.   To do so, we consider a pendulum of length l  =  1m and start the pendulum at a φ0  = 90o   (φ0  = π/2) angle.  We break up the time for 10 full oscillations of the pendulum in roughly 10000 time steps.  Due to numerical errors, the pendulum does not exactly return to φ0  after a full period.  Record and plot the deviation from φ0  at the return point (largest angle) as a function of the number of periods 1, 2, . . . , 10 for both the Euler and the velocity Verlet algorithm.  How does the error increase over time in these two cases?  Try to fit your results for the two cases by a power law.

(b)  (5 pts) Use the velocity Verlet algorithm to obtain the period for one oscillation for l = 1m and initial angles φ0  = π/10, π/4, π/2, π _ 0.1.  Compare with the oscillation period T = 2π ^l/g obtained in the small angle (harmonic) approximation.

(c)  (10 pts) Implement a calculation of the potential energy Epot , the kinetic energy Ekin , and the total energy Etot  = Epot + Ekin .  Record and plot the energies over 10 oscillations of the pendulum using both the Euler and the velocity Verlet algorithm. Compare the results. What do you expect physically?

2.  Motion in a gravitational field                                                                                       (30 points)

Due to gravity, two bodies with masses M and m experience an attractive force

F = _G r

between them.  Here G = 6.674 × 10_ 11  Nm2 /kg2  is the gravitational constant and r the distance between the two bodies.  r is a unit vector connecting the two bodies.  Here we are interested in simulating the motion of bodies due to the gravitational force in two dimensions.

(a)  (20 pts) Write a Python program using Numpy arrays and the velocity Verlet algorithm to

simulate the motion of ten probes in the gravitational field of the sun with mass M = 2 × 1030 kg. The sun is placed at position x = y = 0 and we assume that the probes have such a small mass m that they do not affect the suns position, i.e., the position of the sun is stationary during the simulation. The probes are initially all at the position x = 150 × 109  m and y = 0. We want to understand,  how differences in the initial velocities of the probes affect their trajectories. We therefore consider a distribution of velocities around vx  = 0 and vy  = 30, 000 m/s.   To do so, we add to these velocities random components which can be obtained by using the numpy.random.rand() function.   For example,  the distribution of velocities in x and y direction could have a width of 20, 000 m/s centered around the velocities above, i.e., vx  ∈ [_10000, 10000] and vy  ∈ [20000, 40000]. Plot the trajectories of the probes.