From our previous central force problem, we classified three different types of possible orbits: 1) unbounded orbits, 2) bounded elliptic orbits, and 3) bounded circular orbits. For this assignemnt, we will try to plot some of these orbits by integrating the first integral of motion,
where E is the total energy of the system and the right hand side is simply the sum of the effective potential energy
and the kinetic energy of the system. By rearranging the terms and solving for the radial velocity, we have
Since the angular momentum is also a constant of the motion, we can write
.
This gives us a relation between the differential change dt and the corresponding change dq
.
Substituting this into Eq. (1), we can eliminate dt from the equation. Then, by writing the r dependence explicitly and integrating the equation, we have
,
where r0 and q0 are the initial radial and angular position of the orbit. This equation can then be numerically integated to get the equation of motion for a particular orbit with a given energy E and angular momentum l.
For this problem, we again assume the following form for the interaction potential U(r)
and we take l2/2m=0.01, k=70, a=0.2. Then, do the following:
- Unbounded orbit (E=+25): Taking the following initial condition for the orbit,
numerically integrate Eq. (2) to r=rturning using the Simpson's rule. The orbit will be given in polar coordinate by (r,q) calculated from this procedure. Note that the integrant is singular at r=rturning but the integral with rturning as one of its limit points remains well-defined*. This procedure will give the incoming branch of the orbit and the outgoing branch can be also be easily obtained by symmetry. In particular, the orbit is symmetric with repect to its q coordinate at the turning points. Plot the complete orbit including its incoming and outgoing parts for this case.- Bounded elliptic orbit (E=-25): Similar to the unbounded case, the orbit is invariant under reflection about the apsidal vectors - the vector connecting the center of force and the turning points (see the following figure).
Notice that the two extensions are reflections of the basic orbit about the two blue apsidal vectors. However, for the choice of our U(r), the bounded orbits will not be closed. They will form "spirograph" patterns within an annulus domain defined between rmin and rmax. For this part of the problem, you can first construct the basic orbit by integrating Eq. (2) between rmin and rmax using the Simpson's rule. [Note: Integrals with limits of integration starting or ending at rmin or rmax will be close to improper*. One should consider on how to effectively in dealing with them.] Then, the other extensions can be constructed by reflecting the basic orbits (and other extensions) about the apsidal vectors. For a visually pleasing result, you should plot at least seven or eight extensions.- In the bounded orbit case, will the orbits in general be closed (meaning that will it close on itself after a finite number of loops)?
Your answers are required to be accurate up to single precision.
*This is an example of a treatable improper integral. If one of the limit of the integral is exactly on the turning point, an open formula should be used. However, in our case, the square-root term in the denominator is only approximately zero since the turning point was approximated by another numerical method. Thus, the simple Simpson's rule will work in this case. On the other hand, the integrant still varies greatly near this turning point so that many subdivisions for the Simpson's rule might be needed to get to a desired tolerance. One should note that in plotting the orbits, you have the freedom in choosing the starting points and you can use this freedom to minimize the problems with imporper integrals. We will discuss the proper treatment of singularity in more detail in the next lecture.
Other related links on central force orbits: at Queen's Univ and at UMD
Solution
The implementation of the solution in C can be found here. Since the Simpson's rule will typically require many subdivision if one of the intergal's limits is near a singluar point, an efficient strategy is to start the plotting away from the turning points and move toward them. In this way, only the a few integrals will be close to "improper" (with limits pts close to the turning points). For the unbounded orbit, one may start with r=1 and move toward rturning. The outgoing orbit is symmetric to the incoming orbit with respect to the turning point. In particular, If our incoming orbit is given by the sequence of points {ri,qi} for i=0,...,N, then the outgoing orbit will be given by the sequence of points {rN-i,qN+(qN- qN-i)}. For the bounded case, one should start at a midpoint between rmin and rmax and use two separate Simpson loops to move the orbit toward the respective turning points. If we start plotting from the turning points, all integrals will be close to imporper and the Simpson's Rule might take many steps.
The followings are the graphs for the unbounded and bounded orbits. There are 1000 data points for the basic orbit in each case. The unbounded orbit is with one extension and the bounded orbit is with 20 extensions.
The red circles (ovals) indicate the radial locations of the turning points. The unbounded orbit travels from r=1 toward the center of force and then it scatters away from the center of force at its closest approach at r=rturning (the small red circle). The bounded orbit is forever restricted to move between rmin and rmax (the two red circles). The resulting pattern is similar to a spriograph. However, the bounded orbit is not closed. It will never retrace its own path. This is a result of the Bertrand's theorem which stats that the only types of central force capable of producing closed orbits independent of initial conditions are either an inverse-square force law or a hook-type force law.
As we will learn later, dealing with the differenial equations directly and solving them with a ODE solver will be a more efficient method to plot these orbits.
