MSA-Error Scaling for 2nd-Order Schemes-Ruby
From Micasim
| Project | ACS |
|---|---|
| Project Lead | Piet Hut |
| Project Homepage | ACS |
Contents |
Error Scaling for 2nd-Order Schemes (Ruby Version)
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)
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 10-15 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.
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)
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
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.
