In the previous chapters we saw how slow the first-order forward Euler scheme converged. The leapfrog scheme, being second order, is expected to fare better. Indeed, here are the results, in the form of energy output, output of the last couple positions and velocities, and Figs. 4.1, 4.2.

|gravity> leapfrog1 > leapfrog1_0.01.out Please provide a value for the time step 0.01 Initial total energy E_in = -0.875 Final total energy E_out = -0.87497 absolute energy error: E_out - E_in = 2.98294e-05 relative energy error: (E_out - E_in) / E_in = -3.40907e-05 |gravity> tail -2 !$ tail -2 leapfrog_0.01.out 0.583527 -0.387366 0 1.03799 0.167802 0 0.593822 -0.385632 0 1.02114 0.178871 0 |gravity> leapfrog1 > leapfrog1_0.001.out Please provide a value for the time step 0.001 Initial total energy E_in = -0.875 Final total energy E_out = -0.875 absolute energy error: E_out - E_in = 3.17517e-07 relative energy error: (E_out - E_in) / E_in = -3.62876e-07 |gravity> tail -2 !$ tail -2 leapfrog1_0.001.out 0.590112 -0.362786 0 1.04675 0.203781 0 0.600491 -0.360694 0 1.02914 0.214483 0 |gravity>

We can see how the orbit still changes after each revolution, due to
numerical errors, in Fig. 4.1. However, with a time
step of `0.001`, Fig. 4.2 shows no discernible
errors any more. Unlike the case with the forward Euler integrator,
the last couple output values for position and velocity differ by only
a few percent at most. And finally, it is clear that the leapfrog is
indeed second-order: making the time step ten times smaller reduces
the errors by a factor a hundred. The smaller the time step, the more
accurate this factor of one hundred becomes:

|gravity> leapfrog1 > /dev/null Please provide a value for the time step 0.0001 Initial total energy E_in = -0.875 Final total energy E_out = -0.875 absolute energy error: E_out - E_in = 3.19364e-09 relative energy error: (E_out - E_in) / E_in = -3.64987e-09 |gravity> !! Please provide a value for the time step 0.00001 Initial total energy E_in = -0.875 Final total energy E_out = -0.875 absolute energy error: E_out - E_in = 3.19509e-11 relative energy error: (E_out - E_in) / E_in = -3.65153e-11 |gravity> !! Please provide a value for the time step 0.000001 Initial total energy E_in = -0.875 Final total energy E_out = -0.875 absolute energy error: E_out - E_in = 8.60312e-13 relative energy error: (E_out - E_in) / E_in = -9.83214e-13 |gravity> !! Please provide a value for the time step 0.0000001 Initial total energy E_in = -0.875 Final total energy E_out = -0.875 absolute energy error: E_out - E_in = 1.11888e-12 relative energy error: (E_out - E_in) / E_in = -1.27872e-12

Notice how the scaling first takes this factor 100 to much better than
1% accuracy, but then flattens off when we begin to approach machine
accuracy for double precision (of order ). In the last
case, the cumulative round-off effects are particularly severe. Not
surprising, given the fact that a time step of `dt = 0.000,000,1`
implied a whooping 100,000,000 integration steps, which took a quarter
of an hour to complete. Whereas the scaling from the first two
outputs would have predicted a final error around `3.20e-15` in
absolute energy and around `3.65e-15` in relative energy, the
actual result shows a turn-up of the error size, paradoxically making
the calculation with one hundred million steps less accurate than the
calculation with ten million steps.

The lesson is to always beware of stretching an algorithm beyond its natural area of application. In practice, we will use the leapfrog when modest accuracy is needed, with errors remaining at a level of, say, to or so. If much higher accuracy is needed, it is better to use a fourth-order scheme, as we will see later on in this book.

平成18年10月10日