Problem set #3. Due  .................................

This set revolves around Orbits.

In the original 1987 edition of BT, chapt. 2 (p. 140), the authors proved that the Lagrange points L4 and L5 of the rotating logarithmic potential are ALWAYS stable. This proof was incorrect according to Pfenniger (1990, A\&A 230, 55), and the new edition of BT has been modified.

1. What was wrong with the original proof?

2. Settle any doubts you might still have by numerically solving an initial-value problem of Newton's equation of motion in the corotating frame of reference (BT eq. 3-82) for the appropriately chosen potential (and rotation speed), trying to demonstrate the instability of motion around L4 or L5.

Of course, the initial position of the test particle should not be precisely the L4/L5 point but some small amount away from it, otherwise you may never see much action occur.

The more sophisticated time-integration scheme you use the better. Consider:
Euler's method (very rudimentary, of the type v=v+dt*f, x=x+dt*v, dt=const.),
variable-timestep Euler scheme, where dt for the next step is computed from the requirement that the relative increase in any quantity like v or x has been smaller than a given value TOL, e.g., TOL=0.01),
and a range of Runge- Kutta methods, for instance the 4th order R-K (find it in any numerical analysis textbook, e.g. Numerical Recipes), or the 7th/8th order R-K method (Aarseth and Proszynski).

Try different time-step values and/or different TOLerance parameters, and analyze the accuracy of the results (e.g., the position of the test particle after a certain time, or Jacobi constant conservation). Compare the cpu time needed to achieve the same accuracy with different methods. How much time did you need to implement/debug your method(s)?

This problem explores some issues relevant to numerical simulations using particles. In order to use a grid-based Poisson/force solver, we must:

1. assign the particle density to the mesh,
2. find the potential and numerically find its gradient; OR find directly Fx, Fy, and Fz = force components
3. assign the forces to particles.

• nearest grid point assignement
• linear interpolation (e.g., bi-linear in 2-D)
• cloud-in-cell method with more sophisticated density distributions of a cloud representing the particle. Any finite-support bell-shaped function that covers at most a few cells will do ( Gaussian has infinite support, ie. reaches too far to be practical.)

Analyze analytically and numerically the self-force on one particle with arbitrary position, for instance x=3.1415926, y=x. Is it zero? If not - why? Can you say what is required and/or sufficient for the self-force to vanish?

Then introduce a second particle, for instance at a position with the negative sign of x, y, or both x and y. Compare results with the exact Newtonian force. How does the error depend on the mesh size (dx,dy)?

Compare the accuracy and computational cost of the Poisson (potential) solver vs. that of the Fourier force solver.