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 and , and determine the magnitudes of their separation in space, as well as the magnitude of the difference between their velocities (their separation in velocity space), . The time scale gives an estimate of the minimum time it will take for these two particle to collide. If their relative velocity vector is not closely aligned with their relative position vector and pointing in the same direction, the magnitudes and may not become that much smaller during the next time interval , but certainly the direction of both and will change significantly. In either case, it is important that the integration time step during this period is chosen to be significantly smaller than , 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 ,
with respect to all and 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 and ,
stopping there rather than taking their square roots. The latter
choice would force us to calculate roughly 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
combination that minimizes
will minimize
, 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 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
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
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()`.

平成18年10月10日