4.3 Finding Better Convergence

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>

Figure 4.1: Relative orbit for the first attempt to integrate a two-body system with a leapfrog integrator, with time step $dt = 0.01$

Figure 4.2: Relative orbit for the second attempt to integrate a two-body system with a leapfrog integrator, with time step $dt = 0.001$

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 $10^{-15}$). 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, $10^{-4}$ to $10^{-6}$ 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.