Erica: I think I know what's going on. The original forward Euler algorithm was making terrible errors at pericenter. I'm sure that the backward Euler algorithm would be making similarly terrible errors when the stars pass each other at pericenter. Now the modified Euler scheme works so well because it almost cancels those two errors.

In other words, during pericenter passage, the errors in forward and in backward Euler grow enormously, and the attempt at canceling is relatively less successful. But once we emerge from the close encounter, the attempts at canceling have paid off, and the net result is a much more accurate energy conservation.

Carol: Hmm, that sounds too much like a hand-waving argument to me. I would be more conservative, and just say that second-order algorithms are more complicated to begin with, so I would expect them to have more complex error behavior as well. Your particular explanation may well be right, but can you prove that?

Erica: I'm not sure how to prove it. It is more of a hunch.

Dan: Let's not get too technical here. We want to move stars around, and we don't need to become full-time numerical analysts.

Erica: But I would like to see what happens when we make the time step smaller.

Erica: But the error peaks scale like a second-order algorithm: they have become 100 times less high.

Carol: So the net error after the whole run must have scaled better than second-order. And indeed, look at the output we got when I did the runs: after ten time units, the energy error became a factor thousand smaller, when I decreased the time step by a factor ten!

Dan: Almost too good to be true.

Carol: Well, a second-order scheme is guaranteed to be at least second order; there is no guarantee that it doesn't do better than that. It may be the particular configuration of the orbit that gives us extra error cancellation, for free. Who knows?

Dan: Let's move on and see what the leapfrog algorithm shows us.

Erica: That makes sense, actually. Remember, the leapfrog is time symmetric. Imagine that the energy errors increased in one direction in time. We could then reverse the direction of time after a few orbits, and we would play the tape backward, returning to our original starting point. But that would mean that in the backward direction, the energy errors would decrease. So that would spoil time symmetry: if the errors were to increase in one direction in time, they should increase in the other direction as well. The only solution that is really time symmetric is to let the energy errors remain constant, neither decreasing nor increasing.

Carol: Apart from roundoff.

Dan: But wait a minute, you just argued so eloquently that the leapfrog algorithm should make almost no error, in the long run.

Erica: Yes, either in the long run, or after completing exactly one orbit -- or any integer number of orbits, for that matter. But in our case, we haven't done either. After ten time units, we did not return to the exact same place.

You see, during one orbit, the leapfrog can make all kind of errors. It is only after one full orbit that time symmetry guarantees that the total energy neither increases or decreases. And compared to other algorithms, the leapfrog would do better in the long run, even if you would not stop after an exact integer number of orbits, compared to a non-time-symmetric scheme, such as modified Euler. The latter would just keep building up errors at every orbit, adding the same systematic contribution each time.

Dan: I'd like to see that. Let's integrate for an exact number of orbits.

Erica: Good idea. I also like to test the idea, to see how well it works in practice. But first we have to know how long it takes for our stars to revolve around each other.

Let's see. In the circular case that we have studied before, while we were debugging, we started with the following initial conditions, where I am adding the subscript c for circular:

23.5. Surprisingly High Accuracy

Dan: Let's bring up a calculator on your screen.

Carol: Why not stay with Ruby and use irb? We can use the fact that :

 |gravity> irb
 include Math
 Object
 pi = 2*acos(0)
 3.14159265358979
 t = (4.0/7.0)**1.5 * 2 * pi
 2.7140809410828
 quit

Erica: How about running it ten times longer? I'm curious to see what will happen.

Carol: Let's find another example of an integer number of orbits, close to a multiple of . Here is one: 201 orbits will take a total time , close to . That should be good enough. Here goes:

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.01
 final time = ?
 545.53
   t = 0,   E_kin = 0.125, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 0, (E_tot - E_init) / E_init = -0
   t = 546,   E_kin = 0.126, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 2.21e-08, (E_tot - E_init) / E_init = -2.53e-08

Dan: So now you're starting to bet as well. I think you're right, but let's check:

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.01
 final time = ?
 545.53
   t = 0,   E_kin = 0.125, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 0, (E_tot - E_init) / E_init = -0
   t = 546,   E_kin = 0.0293, E_pot = -0.0886; E_tot = -0.0594
             E_tot - E_init = 0.816, (E_tot - E_init) / E_init = -0.932

Carol: At least for periodic orbits, such as this one. But I must say, I'm impressed too.

Erica: I had wondered why so many people use the leapfrog algorithm. I'm beginning to see that it has some real advantages!

Carol: So what's next. Do we want to start playing with third-order or fourth-order integration schemes?

Dan: I'd rather go beyond the two-body problem, to the real N-body problem. It's time to do start simulating a real star cluster, rather than just one double star!

Previous

ToC

Up

Next