Previous ToC Up Next

## 21.1. A Matter of Time

Carol: So now we have seen how we are making a really big error, when we use forward Euler with a time step that is really too large for a first-order integration scheme. But that's not what we are really interested in. We want to study the behavior of errors in the case where the time steps are small enough to give reasonable orbit pictures.

Erica: Yes, and then we want to compare first-order and second-order integration schemes, to check whether the energy errors scale in different ways. First we should continue to look at forward Euler, but with smaller time steps.

Dan: You know, I really get tired of writing a whole new file, each time we change a parameter, like the size of a time step. Can't we let the program ask for the time step size, so that we can type it in while the program is running?

Carol: Good idea. That is a much cleaner approach. And while we're at it, why not ask for the total duration of the integration too. And that reminds me, I really wasn't very happy with the way we have forced the code to give sparse output, at every interval of \Delta t = 0.01.

Dan: Remind me, what did we do there?

Dan: Ah, yes, of course, how can I forget such a name!

Carol: Well, as the name says, that code took 100,000 steps, during a total time of 10 time units. With output intervals of length 0.01, this means that we needed only 1,000 outputs. In other words, we needed only to print the results once every 100 steps. We did this with the following rather clumsy trick, using the remainder operator %:

```   if i%100 == 99
print(x, "  ", y, "  ", z, "  ")
print(vx, "  ", vy, "  ", vz, "\n")
end
```

I suggest that instead we introduce the time explicitly. Notice that so far we have used a variable dt to keep the value of the time step size, but we have never kept track of the time itself. Let us introduce a variable t to indicate the time that has passed since the start of the integration. We can then specify a desired interval between outputs, dt_out for short. And to keep track of when the next output is scheduled to happen, we can use a variable t_out. Whenever the time t reaches t_out or goes past it, we need to do another output.

Of course, our diagnostics method should now print the value of the time as well. What else do we need to change. The main loop now becomes a test to see whether the time t has passed to or beyond the final time t_end, specified by the user. And after each output statement, t_out has to be incremented to the next scheduled output time, by adding dt_out to its value. Well, that must be about it, yes?

## 21.2. A New Control Structure

Let me open a file euler_energy_try4.rb and type it all in:

``` require "vector.rb"
include Math

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
r2 = r*r
r3 = r2 * sqrt(r2)
a = -r/r3
r += v*dt
v += 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)
```

Dan: What is gets?

Carol: That is a Ruby input statement, short for `get string'. It reads the next line from the command line. So if you type a single value, and then hit the enter key, gets gobbles up the number you have typed, but packaged as a string. For example, when the code asks you for a time step, and you type

``` 0.01
```
then gets returns the string "0.01", made up of four characters. What we really want is a number, and the method to_f is the built-in Ruby way to convert a string into a floating point number; it is short for `(convert) to float'.

## 21.3. Overshooting

Erica: Let's give it the same values as before, to see whether we get the same output.

Carol: This is what we found before:

``` |gravity> ruby euler_energy_try3.rb > euler_energy_try3.out
E_kin = 0.125, E_pot = -1; E_tot = -0.875
E_tot - E_init = 0, (E_tot - E_init) / E_init = -0
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
```
And let us see what our new program gives:

``` |gravity> ruby euler_energy_try4.rb > euler_energy_try4.out
time step = ?
0.01
final time = ?
10
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 = 10,   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
```
Dan: At least the diagnostics output is the same. How about the output files?

Carol: I'll do a diff:

``` |gravity> diff euler_energy_try3.out euler_energy_try4.out
1001a1002
> 7.7019  -6.2835  0  0.81213  -0.57414  0  -1.45027103130336
```
He, that is strange. Our friend diff claims that the two files are identical except for the fact that our latest code produces one more line of output! Let me check it with a word count:

``` |gravity> wc euler_energy_try.out
1001    7007   59973 euler_energy_try3.out
1002    7014   60033 euler_energy_try4.out
2003   14021  120006 ���
```
And what do you know, yes, 1001 lines in euler_energy_try3.out as it should be, moving from times 0 till 10 with steps of 0.01, but why does euler_energy_try4.out have 1002 lines??

Erica: That last program must be taking one more step, beyond time 10. Can you show the last few lines for both output files?

Carol: Sure:

``` |gravity> tail -3 euler_energy_try3.out
7.6775  -6.2662  0  0.81236  -0.57433  0  -1.45027135833743
7.6856  -6.272  0  0.81229  -0.57426  0  -1.45027124898271
7.6937  -6.2777  0  0.81221  -0.5742  0  -1.45027113997184
|gravity> tail -3 euler_energy_try4.out
7.6856  -6.272  0  0.81229  -0.57426  0  -1.45027124898271
7.6937  -6.2777  0  0.81221  -0.5742  0  -1.45027113997184
7.7019  -6.2835  0  0.81213  -0.57414  0  -1.45027103130336
```
Just like diff told us, the last few lines are identical, except for the fact that euler_energy_try4.out shows one extra step. You must be right: it looks like the code didn't know how to stop in time.

## 21.4. Knowing When To Stop

Erica: I wonder why it overshot.

``` while t < t_end
print "t = ", t, " and t_end = ", t_end, "\n"
```

Here we go:

``` |gravity> ruby euler_energy_try5.rb > euler_energy_try5.out
time step = ?
0.01
final time = ?
10
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 = 10,   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
|gravity> tail euler_energy_try5.out
t = 9.95999999999983 and t_end = 10.0
7.6694  -6.2605  0  0.81244  -0.57439  0  -1.45027146803748
t = 9.96999999999983 and t_end = 10.0
7.6775  -6.2662  0  0.81236  -0.57433  0  -1.45027135833743
t = 9.97999999999983 and t_end = 10.0
7.6856  -6.272  0  0.81229  -0.57426  0  -1.45027124898271
t = 9.98999999999983 and t_end = 10.0
7.6937  -6.2777  0  0.81221  -0.5742  0  -1.45027113997184
t = 9.99999999999983 and t_end = 10.0
7.7019  -6.2835  0  0.81213  -0.57414  0  -1.45027103130336
```
Erica: Aha! The problem is roundoff, that explains everything! The time variable t is a floating point variable, and instead of reaching the exact time 10, after 1,000 steps, it comes ever so close, but not quite at the right point. Therefore, when it runs the loop test, it decides that the time has to be incremented by another time step, and it then overshoots.

Carol: That suggests a simple solution. How about testing whether the time has reached not exactly the end time, but close enough? Close enough could mean half a time step. Let's try! And I'll be bold and call the next file euler_energy.rb, in the hope we now get it right. I will write the loop continuation test like this:

``` while t < t_end - 0.5*dt
```

That should do it:

``` |gravity> ruby euler_energy.rb > euler_energy.out
time step = ?
0.01
final time = ?
10
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 = 10,   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
|gravity> diff euler_energy_try3.out euler_energy.out
```
Dan: And it did it. No differences. Congratulations!

Carol: Let me be double sure:

``` |gravity> tail -1 euler_energy_try3.out
7.6937  -6.2777  0  0.81221  -0.5742  0  -1.45027113997184
|gravity> tail -1 euler_energy.out
7.6937  -6.2777  0  0.81221  -0.5742  0  -1.45027113997184
```
Good. So now we have a new tool, allowing us to change two parameters, without having to change the source code each time. Progress!

## 21.5. Linear Scaling

Dan: Now the point of all this was to check whether the energy errors in forward Euler scale linearly with the time step size. Let's try a few values.

Carol: Sure thing. And now that we can control both the time step and the duration of the integration, let's speed things up a bit, and integrate for only 0.1 time unit. Starting with the previous time step we get:

``` |gravity> ruby euler_energy.rb > euler_energy.out
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.874
E_tot - E_init = 0.000626, (E_tot - E_init) / E_init = -0.000715
```
Making the time step ten times shorter, we find:

``` |gravity> ruby euler_energy.rb > euler_energy.out
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 = 6.26e-05, (E_tot - E_init) / E_init = -7.16e-05
```
And making it yet ten times shorter gives:

``` |gravity> ruby euler_energy.rb > euler_energy.out
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 = 6.26e-06, (E_tot - E_init) / E_init = -7.16e-06
```
Dan: Pretty linear, all right.

Carol: Let's jump to a hundred times smaller time step, to see whether the error still becomes a hundred times smaller:

``` |gravity> ruby euler_energy.rb > euler_energy.out
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.26e-08, (E_tot - E_init) / E_init = -7.16e-08
```
Dan: And so it does.

## 21.6. Picture Time

Erica: I'd like to see how the energy error behaves in time, over a few full orbits, and with better accuracy.

``` |gravity> ruby euler_energy.rb > euler_energy.out
time step = ?
0.0001
final time = ?
10
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 = 10,   E_kin = 1.27, E_pot = -2.07; E_tot = -0.8
E_tot - E_init = 0.0749, (E_tot - E_init) / E_init = -0.0856

```
Erica: Ah, yes, that time step was just about short enough to begin to see the intended orbit, without too much drift.

Erica: Even though the orbit behaves a lot better now, it is clear that the energy errors are still being generated mostly around pericenter passage. In between those close encounters, the energy is very well conserved. But whenever the two stars swing around each other, the energy drifts in a systematic and cumulative way.

Dan: Yes, dramatically so! I can see why people like to use individual time steps. If you use thrown in a few more time steps during close encounters, you can get very much more accuracy as a return for investing very little extra computer time.
 Previous ToC Up Next