Previous | ToC | Up | Next |

** Dan**: I am glad that we now got a few second-order integration schemes.
They surely are a lot more efficient than the first-order forward Euler
scheme we started with!

** Erica**: Definitely. For the same amount of computer time, the
accuracy of second-order schemes is much higher. You know, it would
be nice to quantify that notion, to show exactly how accurate each
scheme really is.

** Carol**: We have done something like that already, by checking how the
endpoint of an orbit converged to a specific value, for smaller and
smaller time steps.

** Erica**: Yes, but in that case we always needed two different choices of
the time step, for two different integrations, so that we could compare
the distance between the two end points. I would prefer to use a measure
that tells us how good a single orbit calculation is. And this is indeed
what astronomers do when they compute orbits: they pick a physical quantity
that should be conversed, and they use that to get an impression of
the size of the numerical errors.

** Dan**: What sort of quantities do you have in mind?

** Erica**: The typical conserved quantities for a system of interacting
particles are energy and angular momentum. Of these, energy is a scalar
and angular momentum is a vector. Therefore, for simplicity, people
generally like to measure the change in energy, in order to get an
idea of the errors introduced during orbit integration.

** Dan**: Okay, let's write a method to check energy conservation.
For a system of two particles, how do you write down the total energy?

** Erica**: There are two contributions. There is the energy of motion,
also called kinetic energy. This energy depends only on the speed of
each particle. For a particle with mass and velocity
, the kinetic energy is

** Dan**: Why is there a minus sign?

** Erica**: The gravitational potential energy is normally chosen to be
zero when the two particles are very far away from each other, which
makes sense, since in that case there is almost no gravitational
interaction. Indeed, in our expression above, for
you can see that
.

It is clear from the definition of that the kinetic energy for each particle is always positive or zero. This implies, because the total energy is conserved, that the potential energy has to be zero or negative.

For example, if you place two particles at rest at a very large distance, the kinetic energy is zero and the potential energy is almost zero as well. Then, when the particles start falling toward each other, the kinetic energy gets larger and larger, and therefore more and more positive. The only way that the total energy can be conserved is for the potential energy to become more and more negative.

** Carol**: In our computer programs we have used the one-body representation,
using the relative separation and relative velocity
, rather than the individual positions
and velocities of the particles. So we have to rewrite
your expressions.

** Erica**: Yes, we have to transform the kinetic and potential energies
from the two-body representation to the one-body representation. Let
us start with the kinetic energy. Using Eq. (93,
we get:

In our case, we have decided to use units in which , so the last two expressions simplify to:

** Carol**: That's a bit annoying, to see that factor
coming in. So far, we did not have to specify the masses of the stars.
Our equation of motion for Newtonian gravity, Eq. (40)
contained only the sum of the masses. So when we choose that sum to
be unity, the equation became simply Eq. (42),
and we had gotten rid of any mention of masses.

In other words, whether we had equal masses, , whether we took one mass to be three times as large as the other, , in both cases the orbits would be exactly the same. However, you are now telling us that the total energy will be different for those two cases. In the first case, the factor whereas in the second case, it becomes , a smaller value.

It seems that when we want to measure energy conservation, we have to make an extra choice, for example by specifying the ratio of the two masses.

** Dan**: But the only thing we care about is whether the energy is conserved.
All we want to know whether the term between parentheses in
Eq. (101) remains constant, or almost constant.
Who cares about the funny factor in front?

** Erica**: Well, yes, I basically agree, but let us try to be a bit more
precise, in describing what we do. To make things clear, let me go back
to our earlier description, which still contains the total mass and the
gravitational constant. If you look at the text books, you will find
that they introduce the so-called reduced mass ,
defined as:

** Dan**: Well, name or no name, let's see whether it is actually
conserved reasonably well in our calculations.

** Carol**: Let us start with our simplest vector code,
euler_vector.rb, and let us add some nice diagnostics.
We definitely want to check to what extent the total energy is
constant, but I would also like to see how the kinetic and potential
energy are varying.

To begin with, let me define a method ` energies` which returns
all three energies, kinetic, potential and total, in one array:

def energies(r,v) ekin = 0.5*v*v epot = -1/sqrt(r*r) [ekin, epot, ekin+epot] end

** Erica**: Specific energies, that is.

** Carol**: Yes, but I don't feel that it is necessary to add that to
the method name.

** Dan**: I knew you would get tired of long names!

** Carol**: Only if there is no danger for confusion. In this case,
we're not going to mix specific and absolute energies -- and if we
ever do, we can always make the name longer again.

Next I would like to write a method that prints diagnostics, let me
just call it `print_diagnostics`, that will print out all
three energies.

** Erica**: But what we really need is to know the * deviation* from the original
energy value, to test energy conservation.

** Carol**: You're right. So that means that we had better measure the energy
right from the start, and then remember that value. Let us call the initial
energy `e0`.

** Dan**: That's almost too short a name, for my taste!

** Carol**: I was just trying to see how far I could push your taste!
Now we can use the method above, picking out the last array element
using the ` Array` method ` last`, to find the initial total energy

e0 = energies(r,v).last

before we enter the integration loop.

Now let me think a bit carefully about the layout that
`print_diagnostics` should follow. We probably only need
three significant digits for the relative energy change. That
will be good enough to see how good or bad our energy conservation
is. In Ruby you can use a C like notation to fix the output format,
using expressions like `printf("%.3g", x)` to print the value
of a floating point variable ` x` with three significant digits.

Actually, what I will do is use ` sprintf` instead of ` printf`, which
prints the same information onto a string, rather than directly onto
the output channel. That way, I can use the Ruby ` print` command,
which takes multiple arguments. If ` x` contains the number
, say, then writing

print "x = ", sprintf("%.3g\n", x)gets translated into

print "x = ", "3.14\n"and since

print "x = 3.14\n"so you should see

x = 3.14on the screen, where I have added the new line character

Hmmm, this ought to do it. Here is the whole code:

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(r,v) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} print "\n" end def print_diagnostics(r,v,e0) ekin, epot, etot = energies(r,v) 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 dt = 0.01 e0 = energies(r,v).last print_pos_vel(r,v) print_diagnostics(r,v,e0) 1000.times{ r2 = r*r r3 = r2 * sqrt(r2) a = -r/r3 r += v*dt v += a*dt print_pos_vel(r,v) } print_diagnostics(r,v,e0)

** Dan**: What does ` STDERR` mean, in front of the print statements?

** Carol**: That means that the information will be printed on the standard
error stream. By default, information will be printed on the standard
out stream. These are the two main output streams in Unix.

** Dan**: Why is it called *error stream*?

** Carol**: If you have a lot of output, you want to redirect that to a file.
But if something suddenly goes wrong, you would like to be warned about
that on the screen. You certainly don't want a warning message or error
message to be mixed in with all the other stuff in the output file; you
might never notice it there.

In our case, I would like to use this error channel to report on the behavior of the energies. In fact, we want to determine the energy errors, so it is somewhat appropriate to use the error stream, even though the name suggests that it is normally used to report real errors. But why not? We will use it here to report on small numerical errors.

** Dan**: So you report the values of the various energy contributions only
at the beginning and at the end of the run.

** Carol**: For now, that is good enough. At least it's a start. But let
me check to see whether all this works. We don't need the positions and
velocities for now, so I will redirect those to our waste basket
`/dev/null`

|gravity> ruby euler_energy_try1.rb > /dev/null 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

** Erica**: As it should be: remember, this was our very first run, when we
used a time step that was too big to integrate an elliptic orbit! We made
a huge error at pericenter. In fact, we can now see that the energy changed
sign. We started with a bound orbit, with a total energy that was negative.
But at the end of the integration the energy has become positive!
That means that the particles can escape to infinity.

** Carol**: Why is that?

** Erica**: When the particles are very far away from each other, the potential
energy becomes negligible, and the energy is dominated by the kinetic energy.
Since kinetic energy cannot be negative, such a wide separation is impossible
if the total energy is negative. But for zero or positive total energy, there
is nothing that can prevent the two particles to fly away from each other
completely. And clearly, that is what happened here, as a result of numerical
errors.

** Dan**: Before drawing too many conclusions, we'd better check whether we
still are talking about the same orbit as we did before.

** Carol**: My pleasure. Here is what the old code gave:

|gravity> ruby euler_vector.rb > euler_vector.out |gravity> head -1 euler_vector.out 1 0 0 0 0.5 0 |gravity> tail -1 euler_vector.out 7.6937453936572 -6.27772005661599 0.0 0.812206830641815 -0.574200201239989 0.0And here is what the diagnostics produces, also at the very beginning and end of the output file:

|gravity> ruby euler_energy_try1.rb > euler_energy_try1.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 |gravity> head -1 euler_energy_try1.out 1 0 0 0 0.5 0 |gravity> tail -1 euler_energy_try1.out 7.6937453936572 -6.27772005661599 0.0 0.812206830641815 -0.574200201239989 0.0

** Erica**: We know that our first orbit integration produced large errors,
and we have quantified that by looking at the final energy error, at
the end of our orbit integration. But it would be a lot more instructive
to see how the energy error is growing in time.

** Carol**: Easily done: in file euler_energy_try2.rb, I will modify
our `print_pos_vel` method to include the total energy value as well,
calling it `print_pos_vel_energy` instead:

def print_pos_vel_energy(r,v,e0) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} etot = energies(r,v).last print (etot-e0)/e0 print "\n" end

As you can see, I am printing the energy last, after the positions and velocities. And of course, in the code I'm replacing the old name by the new name in the two invocations, just before entering the loop and at the end of the loop. Let's run it:

|gravity> ruby euler_energy_try2.rb > euler_energy_try2.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 |gravity> head -1 euler_energy_try2.out 1 0 0 0 0.5 0 -0.0 |gravity> tail -1 euler_energy_try2.out 7.6937453936572 -6.27772005661599 0.0 0.812206830641815 -0.574200201239989 0.0 -1.45027113997184

** Carol**: That's easy to fix. In file euler_energy_try3.rb
I will change the frequent output method as follows:

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

This should look better now:

|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 |gravity> head -1 euler_energy_try3.out 1 0 0 0 0.5 0 -0.0 |gravity> tail -1 euler_energy_try3.out 7.6937 -6.2777 0 0.81221 -0.5742 0 -1.45027113997184

** Carol**: Let's figure out what magic incantations ` gnuplot` wants us to
give it, to plot the energy as a function of time. To start with, let
me remember how we have plotted the orbit. Ah yes, we have been using
the fact that gnuplot use the first two columns, if you don't specify
anything otherwise. Instead of relaying on that default choice, let's
plot the orbit again, this time giving gnuplot explicit instructions to
use the data from columns 1 and 2:

|gravity> gnuplot gnuplot> set size ratio -1 gnuplot> plot "euler_energy_try3.out" using 1:2 gnuplot> quitSo this gives figure 43, and indeed, it looks just like before, when we produced figure 18.

Now to make Carol happy, we will plot the values of the total energy, which reside in column 7.

** Dan**: But wait, that is different. First we were plotting ` y` as a function
of ` x`. Now you are going to plot the energy ` E` as a function of what?
Of time, I guess.

** Carol**: Yes, that would be the most obvious choice. And because we are
using constant time steps, that boils down to plotting ` E` as a function
of output number, if we number the output lines successively. And indeed,
gnuplot does have a way to use the output line number as the thing to plot
along the horizontal axis: if you specify the value 0 as a column number,
the output line number will be used.

** Dan**: Ah, that makes sense, and that is easy to remember. If you have
an output line that reads, say, in the first three columns:

8 20 3 4 5 6 9 2 1 ...then it is as if gnuplot itself adds the line numbers to the left:

1 8 20 3 2 4 5 6 3 9 2 1 ...and now the column numbering starts at 0 instead of at 1.

** Carol**: Yes, come to think of it, that must be the reason they
introduced that notation. Well, let me try:

|gravity> gnuplot gnuplot> set size ratio -1 gnuplot> plot "euler_energy_try3.out" using 0:7 gnuplot> quit

** Carol**: Why would that be?

** Erica**: When the two stars come closer, the orbit becomes more curved,
and in addition, the speed becomes larger. So for both reasons, there
is more of a change in the orbit during a constant interval in time.
It would be great if we could use a smaller time step during pericenter
passage, and I'm sure we'll get to that point, later on. But for now,
as long as we are using constant time steps, a higher speed means that
each step will cover a larger distance in space. So we are in a
situation that we are actually taking longer steps in space exactly
there where the orbit is curved most.

** Dan**: Not a good thing to do.

** Erica**: I agree, but it's the simplest thing to do. We can later try
to be more clever.

Previous | ToC | Up | Next |