7.4 Variable Time Steps: a Simple Collision Criterium

There is still one major shortcoming in hermite4.C, namely the use of a constant time step, which is determined at the beginning of the program, and then used unchanged during the whole orbit integration. For a specific application, where we want to follow the orbits of a small number of particles for a relatively short time, this limitation may not seem so terrible. We can simply rerun the orbit calculations for shorter and shorter time steps, until we see that the orbits no longer change significantly. In fact, this is what we have done with every application in part I. In practice, however, this way of using an integration code is very unsatisfactory.

It would be far better to use variable time steps. If none of the particles is particularly close to any of its neighbors, there is no reason to waste computer time on letting the whole system crawl forwards with tiny time steps, just to prepare for a future time where two or more particles do swing by each other in a short time interval. On the other hand, during such brief periods of fast encounter activity, even a very small initial time step may prove to be not small enough: there is no way of knowing in advance how much we have to slow down in order to guarantee good performance of our integrator. It would be far better to leave the decision concerning the size of the time step up to the computer. An autonomous choice of time step allows us to let the computer time migrate to the `hot spots' in time where sudden fast changes demand higher time resolution.

To add a simple form of this type of cleverness is surprisingly simple. Take each pair of particles $i$ and $j$, and determine the magnitudes of their separation in space, $\vert r_{ij}\vert$ as well as the magnitude of the difference between their velocities (their separation in velocity space), $\vert v_{ij}\vert$. The time scale $t_{ij} = \vert r_{ij}\vert/\vert v_{ij}\vert$ gives an estimate of the minimum time it will take for these two particle to collide. If their relative velocity vector ${\bf v}_{ij}$ is not closely aligned with their relative position vector ${\bf r}_{ij}$ and pointing in the same direction, the magnitudes $\vert r_{ij}\vert$ and $\vert v_{ij}\vert$ may not become that much smaller during the next time interval $t_{ij}$, but certainly the direction of both ${\bf r}_{ij}$ and ${\bf v}_{ij}$ will change significantly. In either case, it is important that the integration time step during this period is chosen to be significantly smaller than $t_{ij}$, so that the change per time step in relative position and velocity, both with respect to magnitude and with respect to direction, remain small. That way, our fourth-order integrator will be able to operate in a regime where further shrinking of the time step is guaranteed to give a shrinking of errors that is proportional to the fourth power of the time step size.

We can implement this idea by taking the minimal value of $t_{ij}$, with respect to all $i$ and $j$ combinations. It is not expensive to compute this minimum, since much of the work is already done anyway in the inner loop that computes accelerations and jerks. Computing this minimum in the function acc_and_jerk() as a side effect, we can pass its value back to the main() program. The code hermite5.C does this by adding one extra variable to the list of arguments of acc_and_jerk():

Note that the extra argument is declared as double & coll_time. Here the symbol & indicates that the variable coll_time will be passed to the fuction acc_and_jerk() by reference, not by value. This means that we can change coll_time inside the function, while the change persists after we leave the function and return to the main() function that called the function acc_and_jerk(). The variable name coll_time_sq, declared in the first line of acc_and_jerk(), is short hand for ``collision time square''. It is much cheaper and natural to calculate the squares of the magnitudes $\vert r_{ij}\vert$ and $\vert v_{ij}\vert$, stopping there rather than taking their square roots. The latter choice would force us to calculate roughly $N^2$ square roots, which are expensive to calculate, since they take much more time than additions or multiplications on most machines. All we want to know anyway is the minimum of their ratios. The same $\{i, j\}$ combination that minimizes $\vert r_{ij}\vert/\vert v_{ij}\vert$ will minimize $\vert r_{ij}\vert^2/\vert v_{ij}\vert^2$, and once we have determined that minimum, we perform the square root operation at the end of the function acc_and_jerk().

Before we enter the $\{i, j\}$ loop, we have to set coll_time_sq to an arbitrary value that is high enough to guarantee that it exceeds the minimum value that we want to determine. Setting it to $10^{300}$ is surely overkill, but better safe than sorry, and according to the IEEE standard for double precision floating point numbers, this value does not yet lead to overflow. Inside the loop, we let coll_time_sq become shorter and shorter, through the assignment:

coll_time_sq = coll_est_sq;

each time when we find that the estimated collision time for the $\{i, j\}$ pair is shorter than the best estimate so far for coll_time_sq. The only other modification to acc_and_jerk() concerns the computation of v2, the square of the distance in velocity space between particles i and j.

The second function, energy(), remains unchanged. In the main() part of the program, the changes are minor, as can be seen from the result of doing a diff on the main() functions of hermite4.C and hermite5.C:

61c71
<     double dt = atof(argv[1]);
---
>     double dt_param = atof(argv[1]);
81c91,93
<     acc_and_jerk(m, r, v, a, jk, n);
---
>     double coll_time;
>     acc_and_jerk(m, r, v, a, jk, n, coll_time);
>     double dt = dt_param * coll_time;
102c114
<         acc_and_jerk(m, r, v, a, jk, n);
---
>         acc_and_jerk(m, r, v, a, jk, n, coll_time);
110a123,124
> 	dt = dt_param * coll_time;

Instead of reading in the pregiven value for dt from the first argument on the command line, we now read in the value of dt_param, the factor with which we will multiply the estimate for the smallest collision time, in order to determine the next time step, as can be seen in the middle as well as in the very last line of changes reported by diff. The only other changes are the declaration of double coll_time, and the addition of that variable to the list of arguments in the two calls of the function acc_and_jerk().