19. Time Reversibility

19.1. Long Time Behavior

Let me take the modified Euler code, copying it from euler_modified_vector.rb to euler_modified_long_time.rb. I will let the code run ten times longer, by changing the loop defining line to:

 10000.times{                           

Carol: I'm glad you're getting the hang of using long names. Thank you!

Dan: My pleasure. But see, I did still abbreviate a bit: I could have left the word vector in, but that really would have made the name too long, for my taste.

On a more important topic, I really don't like having different files lying around that are almost the same, except for just one extra 0 in one line.

Carol: We'll have to do something about that. I had already been thinking about introducing command line arguments.

Erica: What does that mean?

Carol: We really would like to specify the number of steps on the command line, as an argument. It would be much better if we could take the program euler_modified_vector.rb and run it for 10,000 steps, simply by invoking it as

 |gravity> ruby euler_modified_vector.rb -n 10000

 |gravity> ruby euler_modified_vector.rb -t 100

Dan: I would like that much better! Let's put that on our to-do list. But for now, let me finish my long time behavior test. I'll write leapfrog_long_time.rb, by modifying leapfrog.rb in the same way, to take 10,000 steps:

 10000.times{                  

Our expectation would be that modified Euler will completely screw up, while the leapfrog will keep behaving relatively well. Let's see what will happen!

First I will make a picture of the long run for modified Euler, in figure 39:

 |gravity> ruby euler_modified_long_time.rb > euler_modified_long_time.out

 |gravity> ruby leapfrog_long_time.rb > leapfrog_long_time.out

19.3. Testing Time Symmetry

I will call the new code leapfrog_backward.rb, which is the same as the old code leapfrog_long_time.rb, except that I have replaced the original loop by the following two loops:

 10000.times{
   v += 0.5*a*dt
   r += v*dt
   a = acc(r)
   v += 0.5*a*dt
 }
 
 dt = -dt
 10000.times{
   v += 0.5*a*dt
   r += v*dt
   a = acc(r)
   v += 0.5*a*dt
   print_pos_vel(r,v)
 }

I will plot the backward trajectory in figure 41:

 |gravity> ruby leapfrog_backward.rb > leapfrog_backward.out

Erica: Ah, yes, but that's precisely the point. The stars are retracing their steps so accurately, we can't see the difference!

Dan: Let's check how close the stars reach their point of departure, after their long travel:

 |gravity> ruby leapfrog_backward.rb | tail -1
 0.999999999999975  -1.06571782648723e-12  0.0  2.12594872261995e-12  0.500000000000013  0.0  

 r = [1, 0, 0].to_v      
 v = [0, 0.5, 0].to_v    

19.4. Two Ways to Go Backward

Dan: No, I've gone back in time, using negative time steps, while leaving everything else the same, including the sign of the velocity. I could instead have reversed the direction of the velocity, while leaving the time step the same. That would mean that time would keep flowing forward, but the stars would move in the opposite direction, from time to . Let me try that too, why not, in leapfrog_onward.rb. This code is the same as leapfrog_backward.rb, with the only difference being the one line in between the two loops, which now reads:

 v = -v                  

In this case, the final position and velocity are:

 |gravity> ruby leapfrog_onward.rb | tail -1
 0.999999999999975  -1.06571782648723e-12  0.0  -2.12594872261995e-12  -0.500000000000013  -0.0  

Erica: It is remarkable how close we come to the starting point. And yet, it is not exactly the starting point.

Carol: The small deviations must be related to roundoff. While the algorithm itself is strictly time symmetric, in the real world of computing we typically work with double precision numbers of 64 bits. This means that floating point numbers have a finite precision, and that any calculation in floating point numbers will be rounded off to the floating point number that is closest to the true result. Since the rounding off process is not time symmetric, it will introduce a slight time asymmetry.

19.5. Testing a Lack of Time Symmetry

I will do the same as what we did for the leapfrog. I will call the new code euler_modified_backward.rb, which is the same as the old code euler_modified_long_time.rb, except that I have again replaced the original loop by these two loops:

 10000.times{
   r1, v1 = step_pos_vel(r,v,dt)
   r2, v2 = step_pos_vel(r1,v1,dt)
   r = 0.5 * ( r + r2 )
   v = 0.5 * ( v + v2 )
 }
 
 dt = -dt
 10000.times{
   r1, v1 = step_pos_vel(r,v,dt)
   r2, v2 = step_pos_vel(r1,v1,dt)
   r = 0.5 * ( r + r2 )
   v = 0.5 * ( v + v2 )
   print_pos_vel(r,v)
 }

I will plot the backward trajectory in figure 42:

 |gravity> ruby euler_modified_backward.rb > euler_modified_backward.out