This week's assignment is a very interesting one! As a first step, reading the wiki article on Lagrange points would be a good start. For the mathematically macho amongst you, attempting to derive their location using Newtonian mechanics is a nice exercise. A potential formalism (rather than forces) in a rotating frame of reference is mathematically neater, with a result illustrated in the wiki article. An approach using Lagrangian or Hamiltonian dynamics makes the conclusion pretty straight forward. Such treatments can be found in any mechanics textbook.
One interesting generalisation is to 3D - to permit motions out of the plane (for the massless particle). One interesting result of this the Kozai effect, wherein moons of planets not in the ecliptic have highly unstable orbits.
In the 2D plane case, the n-body problem has a small number of stable solutions, mostly involving orbital resonances, an area of active research and the largest contributing factor to the demotion of Pluto.
This is the restricted 3-body problem mentioned above, wherein one of the bodies has negligible mass. If you derive (for fun) a potential map that looks similar to this, or just borrow it, then you can see the stability of each Lagrange point. In particular, the coaxial points are unstable along the radius. NASA still uses them to place research satellites like SOHO, though in practice it orbits them in a plane perpendicular to the radius. In order to illustrate this orbit and demonstrate its stability, the code will need to be generalised to three dimensions. For well designed code, this is a matter of changing a few numbers.
The second point of interest concerns the Trojan points L4 and L5. As can be seen in the potential diagram, these points are actually unstable in all directions - how then do asteroids accumulate out there? This is something Casey will discuss during the section on Nov 15, but if it interests you, write a sentence in your solution. This is also discussed on the wiki article if you get stumped.
For more general simulations of general 3 body problems (even in 2D), it may be necessary to implement adaptive time steps. How is this done? In general, the next step can be computed using, eg h = 0.1 and h = 0.05, and the results compared. If they are the same, the step is accepted and the time step can be increased. If they are more different than some predefined tolerance, then the step is rejected and the time step reduced, usually by a factor of 2. For the comparison, the previous time steps can be cached to avoid calculating complex things twice (though in this case it doesn't matter too much). For these purposes, it is recommended to write a separate module to manage the adaptive time stepping, and allow it to be turned on or off as required. I'm looking for a degree of creativity in this regard, though of course books have been written on the subject. In general explicit time steppers require smaller step sizes than implicit ones. There is an important measure of the accuracy of a time stepper which gives as a fraction the degree that can be stepped forward for a given spatial resolution, though it is more relevant to solving of field PDEs rather than point-like particles.
Adaptive time stepping is required since the acceleration varies greatly in non circular orbits depending on the phase. Since some parts will accelerate more than others, the code should be designed so that all variables update below the critical error.
Animation is possible in python, though I'm not sure how. Try google. If you have particularly awesome data, animate it in mathematica, export it as an mpeg, add your favourite music and upload it to youtube, then send me the link.
This week's set will test your organisation and development skills, as it requires the construction and testing of three interacting pieces of code.
1) The first is a function or class which tells you where the planets are at any given time. The assignment is quite detailed in how this could be done, but the important thing is that it's reliable, testable, and can output data nicely.
2) The integrator for the trajectory of voyager 2, based on some initial conditions, such as time of launch, velocity of launch, and angle of launch. For purposes of simplicity, you can remove Earth from part 1, or better yet, set its mass to zero.
3) Some sort of optimiser to locate the conditions that are perfect for hitting (or rather, just missing) each of the planets. Casey's section got a quick spiel on MCMC integrators using the Metropolis Hastings Algorithm, a random number based search algorithm that is versatile, functional in high dimensional spaces, extremely simple to write and modify, and can be automated. While it is much less efficient than certain dedicated optimisation algorithms, it can readily be set to run overnight. More details below.
GRAPHING: Graphical presentation of information is becoming more and more important in this class and in science in general as information becomes more complicated and more dense. Ideas to consider and think about when planning a graph include clarity, data/ink ratio, optimal angle, and so on. You can develop your own heuristics for a good graph. For reference, here are some nice examples. If you want to do funky animations and pyplot is no longer cutting it, feel free to use a graphing program of your choice.
Metropolis Hastings Algorithm: Although any level of complexity can be added to the MH Algorithm, its core idea is presented here. We initialise with a set of parameters which define a position in phase space, eg params = (launch velocity, launch angle, launch day). Then we define a cost function that takes params as an argument and delivers a single number to be minimised. For example, distance from Jupiter (or Saturn…), with some high cost given if collision with Jupiter occurs. The algorithm then works as follows.
cost = CostFunction(params)
paramstar = params + random increment (this part can be tricky)
coststar = CostFunction(paramstar)
if cost/coststar>1, then params = paramstar (* see comment)
wash rinse repeat
*Comment: to enable a broader search, replace 1 with randomreal[0,1] and raise the cost ratio to an appropriate power to keep the acceptance ratio at about 25%. Obviously if it accepts every offer, it'll never home in on the optimal position in phase space. Conversely, exclusively downhill searches have the potential to become trapped in local minima.
Additional concerns of automation include building in stages to allow optimisation for each successive planet, often with far smaller random increments of the initial parameters.