Previous | ToC | Up | Next |

** Bob**: Time to run the driver for our new integrator,
which will give us the energy conservation diagnostics
we have been talking about for so long now! Let us start
with the same variables we used for our very first forward
Euler run:

dt = 0.01 # time step dt_dia = 1 # diagnostics printing interval dt_out = 1 # output interval dt_end = 1 # duration of the integration

We concluded that the position was probably right to within a percent or so, from comparing it with a more accurate integration. I'm curious whether the energy check gives us a similar answer.

|gravity> ruby test.rb < euler.in dt = 0.01 dt_dia = 1 dt_out = 1 dt_end = 1 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 1, after 100 steps : E_kin = 0.838 , E_pot = -1.7, E_tot = -0.867 E_tot - E_init = 0.00822 (E_tot - E_init) / E_init =-0.0094 1.0000000000000000e+00 4.4322956434140903e-01 3.8421555868663582e-01 -1.2943653128939152e+00 2.5812000666903590e-02

** Bob**: That was just pure luck, my argument was only an order of
magnitude estimate.

** Alice**: As we argued, a first-order method like the forward Euler
should converge linearly: making the step size ten times smaller
should decrease the error size also by roughly a factor of ten.
Shall we try?

dt = 0.001 # time step dt_dia = 1 # diagnostics printing interval dt_out = 1 # output interval dt_end = 1 # duration of the integration

** Bob**: So we expect an error of now. Let's see.

|gravity> ruby test.rb < euler.in dt = 0.001 dt_dia = 1 dt_out = 1 dt_end = 1 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 1, after 1000 steps : E_kin = 0.864 , E_pot = -1.74, E_tot = -0.874 E_tot - E_init = 0.000825 (E_tot - E_init) / E_init =-0.000943 1.0000000000000000e+00 4.3302178053157708e-01 3.7860453335416983e-01 -1.3147926449876397e+00 7.1843385614138713e-03

dt = 0.00001 # time step dt_dia = 1 # diagnostics printing interval dt_out = 1 # output interval dt_end = 1 # duration of the integration

making the time step smaller by a jump of a factor of a hundred.

|gravity> ruby test.rb < euler.in dt = 1.0e-05 dt_dia = 1 dt_out = 1 dt_end = 1 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 1, after 100000 steps : E_kin = 0.867 , E_pot = -1.74, E_tot = -0.875 E_tot - E_init = 8.25e-06 (E_tot - E_init) / E_init =-9.43e-06 1.0000000000000000e+00 4.3186966438051116e-01 3.7796470836807455e-01 -1.3171480879236508e+00 5.0327840503739613e-03

** Alice**: Let's follow your earlier suggestion to do a longer integration,
for ten time units. And let's use your flexible output options. I would
like to see four energy diagnostics, but at most one final position and
velocity output, so that we don't get too distracted.

** Bob**: That means:

dt = 0.01 # time step dt_dia = 2.5 # diagnostics printing interval dt_out = 10 # output interval dt_end = 10 # duration of the integration

I have given `dt_out` a value larger than the total run time
`dt_end`, so we won't get the regular particle output, only
the diagnostics.

|gravity> ruby test.rb < euler.in dt = 0.01 dt_dia = 2.5 dt_out = 10 dt_end = 10 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 2.5, after 250 steps : E_kin = 0.936 , E_pot = -0.54, E_tot = 0.395 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45 at time t = 5, after 500 steps : E_kin = 0.603 , E_pot = -0.209, E_tot = 0.394 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45 at time t = 7.5, after 750 steps : E_kin = 0.529 , E_pot = -0.135, E_tot = 0.394 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45 at time t = 10, after 1000 steps : E_kin = 0.495 , E_pot = -0.101, E_tot = 0.394 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45 1.0000000000000000e+00 7.6937453936571965e+00 -6.2777200566159870e+00 8.1220683064181454e-01 -5.7420020123998905e-01

** Bob**: Ah yes, this was the run where the system exploded, with the
relative distance increasing to a value of about ten.

** Alice**: And look, the velocity vector is almost parallel to the position
vector, and pointing in the same direction. Clearly, the two bodies are
escaping from each other.

** Bob**: And that explains why the energy error is frozen: there is no longer
much interaction going on between the particles. The numerical explosion
must have happened at pericenter, as we deduced earlier. Let's have a
closer look at what happens in the first part:

dt = 0.01 # time step dt_dia = 0.5 # diagnostics printing interval dt_out = 10 # output interval dt_end = 2.5 # duration of the integration

A four times shorter run, with five times more frequent diagnostics.

|gravity> ruby test.rb < euler.in dt = 0.01 dt_dia = 0.5 dt_out = 10 dt_end = 2.5 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 0.5, after 50 steps : E_kin = 0.231 , E_pot = -1.1, E_tot = -0.872 E_tot - E_init = 0.0033 (E_tot - E_init) / E_init =-0.00377 at time t = 1, after 100 steps : E_kin = 0.838 , E_pot = -1.7, E_tot = -0.867 E_tot - E_init = 0.00822 (E_tot - E_init) / E_init =-0.0094 at time t = 1.5, after 150 steps : E_kin = 3.37 , E_pot = -2.95, E_tot = 0.417 E_tot - E_init = 1.29 (E_tot - E_init) / E_init =-1.48 at time t = 2, after 200 steps : E_kin = 1.26 , E_pot = -0.865, E_tot = 0.398 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45 at time t = 2.5, after 250 steps : E_kin = 0.936 , E_pot = -0.54, E_tot = 0.395 E_tot - E_init = 1.27 (E_tot - E_init) / E_init =-1.45

** Bob**: Let's check:

dt = 0.01 # time step dt_dia = 10 # diagnostics printing interval dt_out = 0.25 # output interval dt_end = 1.5 # duration of the integration

I'll ask only for particle output, once every quarter time unit:

|gravity> ruby test.rb < euler.in dt = 0.01 dt_dia = 10 dt_out = 0.25 dt_end = 1.5 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 1.0000000000000000e+00 9.6984896296142042e-01 1.2383206879933301e-01 -2.5280864394541708e-01 4.8458243922783906e-01 1.0000000000000000e+00 8.7455226567482147e-01 2.3948976867782423e-01 -5.2597914456085060e-01 4.3083603675399151e-01 1.0000000000000000e+00 7.0559656638283363e-01 3.3480149737403997e-01 -8.5114538471899592e-01 3.1158537532636715e-01 1.0000000000000000e+00 4.4322956434140903e-01 3.8421555868663582e-01 -1.2943653128939152e+00 2.5812000666903590e-02 1.0000000000000000e+00 4.3573377318392045e-02 2.9481575851334291e-01 -1.9466715027999648e+00 -1.1087280827004695e+00 1.0000000000000000e+00 -1.1742854070954888e-01 -3.1759379188140036e-01 1.2444634633070424e+00 -2.2783744800811112e+00

** Bob**: Yes, we'll get to graphics soon. But let's do some real lab bench
work here. Ah! You see, at time 1.25 the separation of the particles is
clearly smaller than at other times. That must have been the periastron
passage, roughly. So the period of the orbit must be 2.5, give or take.

** Alice**: Let's check!

** Bob**: We can only check when we use smaller step sizes, since this one
here is numerically exploding. Let's try this:

dt = 0.001 # time step dt_dia = 2.5 # diagnostics printing interval dt_out = 2.5 # output interval dt_end = 2.5 # duration of the integration

A ten times smaller step size, and both types of output at time 2.5:

|gravity> ruby test.rb < euler.in dt = 0.001 dt_dia = 2.5 dt_out = 2.5 dt_end = 2.5 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 2.5, after 2500 steps : E_kin = 0.202 , E_pot = -0.844, E_tot = -0.643 E_tot - E_init = 0.232 (E_tot - E_init) / E_init =-0.266 1.0000000000000000e+00 1.1300052574190800e+00 -3.5392829077199284e-01 5.6745700562449897e-01 2.8574832366093905e-01

1 1 0 0 0.5

dt = 0.0001 # time step dt_dia = 2.5 # diagnostics printing interval dt_out = 2.5 # output interval dt_end = 2.5 # duration of the integration

|gravity> ruby test.rb < euler.in dt = 0.0001 dt_dia = 2.5 dt_out = 2.5 dt_end = 2.5 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 2.5, after 25000 steps : E_kin = 0.144 , E_pot = -0.994, E_tot = -0.85 E_tot - E_init = 0.0255 (E_tot - E_init) / E_init =-0.0291 1.0000000000000000e+00 9.9746583635490227e-01 -1.3247708073972395e-01 2.6197303614764500e-01 4.6897087837512991e-01

dt = 0.00001 # time step dt_dia = 2.5 # diagnostics printing interval dt_out = 2.5 # output interval dt_end = 2.5 # duration of the integration

|gravity> ruby test.rb < euler.in dt = 1.0e-05 dt_dia = 2.5 dt_out = 2.5 dt_end = 2.5 at time t = 0, after 0 steps : E_kin = 0.125 , E_pot = -1, E_tot = -0.875 E_tot - E_init = 0 (E_tot - E_init) / E_init =-0 at time t = 2.5, after 250000 steps : E_kin = 0.143 , E_pot = -1.02, E_tot = -0.872 E_tot - E_init = 0.00257 (E_tot - E_init) / E_init =-0.00294 1.0000000000000000e+00 9.7908881919480084e-01 -1.0885124045716744e-01 2.2087963584449741e-01 4.8637780261985308e-01

Previous | ToC | Up | Next |