Previous ToC Up Next

22. Error Scaling for 2nd-Order Schemes

22.1. Modified Euler

Dan: Let's have a look at our second-order integration schemes. If I understand things correctly, they are supposed to improve energy quadratically, right? When we make the time step ten times smaller, the energy error should become one hundred times smaller.

Carol: That's the expectation, yes. But first I have to write the codes. I will start with the modified Euler algorithm, for which we had written the vector version of the code in file euler_modified_vector.rb. I will open a new file euler_modified_energy.rb, and add the same type of energy output statements and diagnostics as we have done for the forward Euler case. Here we go: 

 require "vector.rb"
 include Math
 
 def step_pos_vel(r,v,dt)
   r2 = r*r
   r3 = r2 * sqrt(r2)
   a = -r/r3
   [r + v*dt, v + a*dt]
 end
 
 def energies(r,v)
   ekin = 0.5*v*v
   epot = -1/sqrt(r*r)
   [ekin, epot, ekin+epot]
 end
 
 def print_pos_vel_energy(r,v,e0)
   r.each{|x| printf("%.5g  ", x)}
   v.each{|x| printf("%.5g  ", x)}
   etot = energies(r,v).last
   print (etot-e0)/e0
   print "\n"
 end
 
 def print_diagnostics(t,r,v,e0)
   ekin, epot, etot = energies(r,v)
   STDERR.print "  t = ", sprintf("%.3g, ", t)
   STDERR.print "  E_kin = ", sprintf("%.3g, ", ekin)
   STDERR.print "E_pot = ", sprintf("%.3g; ", epot)
   STDERR.print "E_tot = ", sprintf("%.3g\n", etot)
   STDERR.print "            E_tot - E_init = ", sprintf("%.3g, ", etot-e0)
   STDERR.print "(E_tot - E_init) / E_init = ", sprintf("%.3g\n", (etot-e0)/e0)
 end
 
 r = [1, 0, 0].to_v
 v = [0, 0.5, 0].to_v
 e0 = energies(r,v).last
 
 t = 0
 t_out = 0
 dt_out = 0.01
 STDERR.print "time step = ?\n"
 dt = gets.to_f
 STDERR.print "final time = ?\n"
 t_end = gets.to_f
 
 print_pos_vel_energy(r,v,e0)
 t_out += dt_out
 print_diagnostics(t,r,v,e0)
 
 while t < t_end - 0.5*dt
   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 )
   t += dt
   if t >= t_out
     print_pos_vel_energy(r,v,e0)
     t_out += dt_out
   end
 end
 print_diagnostics(t,r,v,e0)

22.2. Energy Error Scaling

Dan: That all looks pretty straightforward.

Carol: And now it is just a matter of making the time steps smaller and smaller: 

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.01
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -8.33e-08, (E_tot - E_init) / E_init = 9.52e-08

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -7.19e-10, (E_tot - E_init) / E_init = 8.22e-10

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.0001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -7.07e-12, (E_tot - E_init) / E_init = 8.08e-12

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.00001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -4.25e-14, (E_tot - E_init) / E_init = 4.86e-14

 |gravity> ruby euler_modified_energy.rb > /dev/null
 time step = ?
 0.000001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -1.41e-13, (E_tot - E_init) / E_init = 1.61e-13
Erica: Up till the last run, it looked almost too good to be true. We must have hit roundoff, I guess.

Carol: Well, yes, with double precision you can't get much further than in relative accuracy for a single calculation. I'm surprised we got as close as we did. Most of the roundoff errors must have cancelled, in the 10,000 steps we took in the next to last integration. But in the last run, where we took 100,000 steps, we accumulated more roundoff errors. When you are adding more steps, you'll get more roundoff, no matter how accurate each individual step may be.

Dan: But wait, just one thing: we haven't checked yet whether we are still getting the same results as before.

Carol: Ah, yes, safety first! The old code gave: 

 |gravity> ruby euler_modified_vector.rb | tail -1
 0.400020239524913  0.343214474344616  0.0  -1.48390077762002  -0.0155803976141248  0.0
and we should get the same result for our new code, if we give it the same parameters: 

 |gravity> ruby euler_modified_energy.rb | tail -1
 0.400020239524913  0.343214474344616  0.0  -1.48390077762002  -0.0155803976141248  0.0
Dan: All is well.

22.3. Leapfrog

Carol: Finally, time to let the leapfrog algorithm tell us whether it is really a second-order algorithm as well. I will start with leapfrog.rb. I will open a new file leapfrog_energy.rb, and again I will add the same type of energy output statements and diagnostics. Here it is: 

 require "vector.rb"
 include Math
 
 def acc(r)
   r2 = r*r
   r3 = r2 * sqrt(r2)
   -r/r3
 end
 
 def energies(r,v)
   ekin = 0.5*v*v
   epot = -1/sqrt(r*r)
   [ekin, epot, ekin+epot]
 end
 
 def print_pos_vel_energy(r,v,e0)
   r.each{|x| printf("%.5g  ", x)}
   v.each{|x| printf("%.5g  ", x)}
   etot = energies(r,v).last
   print (etot-e0)/e0
   print "\n"
 end
 
 def print_diagnostics(t,r,v,e0)
   ekin, epot, etot = energies(r,v)
   STDERR.print "  t = ", sprintf("%.3g, ", t)
   STDERR.print "  E_kin = ", sprintf("%.3g, ", ekin)
   STDERR.print "E_pot = ", sprintf("%.3g; ", epot)
   STDERR.print "E_tot = ", sprintf("%.3g\n", etot)
   STDERR.print "            E_tot - E_init = ", sprintf("%.3g, ", etot-e0)
   STDERR.print "(E_tot - E_init) / E_init = ", sprintf("%.3g\n", (etot-e0)/e0)
 end
 
 r = [1, 0, 0].to_v
 v = [0, 0.5, 0].to_v
 e0 = energies(r,v).last
 
 t = 0
 t_out = 0
 dt_out = 0.01
 STDERR.print "time step = ?\n"
 dt = gets.to_f
 STDERR.print "final time = ?\n"
 t_end = gets.to_f
 
 print_pos_vel_energy(r,v,e0)
 t_out += dt_out
 print_diagnostics(t,r,v,e0)
 
 a = acc(r)
 while t < t_end - 0.5*dt
   v += 0.5*a*dt
   r += v*dt
   a = acc(r)
   v += 0.5*a*dt
   t += dt
   if t >= t_out
     print_pos_vel_energy(r,v,e0)
     t_out += dt_out
   end
 end
 print_diagnostics(t,r,v,e0)

22.4. Another Error Scaling Exercise

Carol: I'll just use the same parameters as I did while testing the modified Euler scheme, for making the time steps smaller and smaller: 

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.01
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 1.19e-07, (E_tot - E_init) / E_init = -1.36e-07

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 1.19e-09, (E_tot - E_init) / E_init = -1.36e-09

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.0001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 1.19e-11, (E_tot - E_init) / E_init = -1.36e-11

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.00001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 1.15e-13, (E_tot - E_init) / E_init = -1.31e-13

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.000001
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -6.22e-15, (E_tot - E_init) / E_init = 7.11e-15

22.5. Roundoff Kicks In

Erica: An amazing accuracy, and that after 100,000 steps! What happened? I would have thought that the cumulative effects of 100,000 roundoff errors would have spoiled the fun.

Carol: We were probably just lucky in the way the roundoff errors canceled. Notice that the energy error at first got 100 times smaller, each time we made the time step 10 times smaller, as it should for a second-order algorithm. But in the last round, the improvement was a lot less than a factor 100.

We must be really close to the roundoff barrier now. Let me just make the time step a factor two smaller. That should make the total error grow again. Wanna bet?

Dan: No, but I do wanna see whether you're right.

Carol: Here we go; 

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.0000005
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = 4.66e-15, (E_tot - E_init) / E_init = -5.33e-15
Dan: I should have accepted your challenge, and made a bet against you!

Carol: I must say, I'm surprised that the roundoff errors cancel so well. But this just can't go on. If shrink the time step factor by another factor of two. 

 |gravity> ruby leapfrog_energy.rb > /dev/null
 time step = ?
 0.00000025
 final time = ?
 0.1
   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 = 0.1,   E_kin = 0.129, E_pot = -1; E_tot = -0.875
             E_tot - E_init = -2.74e-14, (E_tot - E_init) / E_init = 3.13e-14
So there! Now the errors are finally accumulating enough to show a worse performance.

Dan: But you didn't feel confident enough to ask us to bet, this time.

Carol: I should have! And yes, before you ask me, let us check whether we still get the same output as before. What we got was: 

 |gravity> ruby leapfrog.rb | tail -1
 0.583527377458303  -0.387366076048216  0.0  1.03799194001953  0.167802127213742  0.0
And what our new code gives is: 

 |gravity> ruby leapfrog_energy.rb | tail -1
 0.583527377458303  -0.387366076048216  0.0  1.03799194001953  0.167802127213742  0.0
Good.
Previous ToC Up Next