Previous | ToC | Up | Next |

** Carol**: This is about the simplest thing we could possibly do,
for the one-body problem: starting with a circular orbit, and then
taking only one small step. Here we go . . .

|gravity> ruby euler_try_circular_step.rb 1 0 0 0 1 0 1.0 0.01 0.0 -0.01 0.99 0.0

** Erica**: Simple, yes, but correct? Let's compute the numbers by
hand, to see whether our program gave the right answers. We start
with

** Dan**: That is encouraging!
Now what about the second half of the second output lines?

** Erica**: To compute the new velocity, we have to first compute the
acceleration vector. We can use Eq. (42), which
I'll copy here once again:

** Dan**: Spot on! We should have gotten , but
we somehow wound up with . So that's were
the bug is, in the * y* component of the new velocity, which should be
, but somehow isn't.

** Erica**: Easy to check: here is where we compute the new value of
, which we call `vy`:

vx += ax*dt vy += ax*dt vz += az*dt

Ooops!! A typo. How silly. No wonder we got the wrong answer for . Let me correct it right away, and write it out as a new file, euler_circular_step.rb:

vx += ax*dt vy += ay*dt vz += az*dt

I'm curious to see whether now everything will be all right:

|gravity> ruby euler_circular_step.rb 1 0 0 0 1 0 1.0 0.01 0.0 -0.01 1.0 0.0

** Carol**: Yes, we have now verified that we got the right result after
one step.

** Dan**: Great! Let's go back to our original code, correct the bug,
and we'll be in business.

** Carol**: Only if the first bug we caught will be the last bug.
Don't be so sure! We may well have made another mistake somewhere
else.

** Dan**: Oh, Carol, you're too pessimistic. I bet everything will be
fine now!

** Erica**: We'll see. Here is the same typo in euler_try.rb:

vx += ax*dt vy += ax*dt vz += az*dt

and here is the corrected program, which I will call euler.rb in the spirit of Dan's optimism:

vx += ax*dt vy += ay*dt vz += az*dt

I'll run the new code:

|gravity> ruby euler.rb > euler.outAnd here is the plot, in fig 18:

|gravity> gnuplot gnuplot> set size ratio -1 gnuplot> plot "euler.out" gnuplot> quit

** Dan**: No comment.

** Erica**: Maybe we should go back to the circular orbit. We tried to
take a single step there, and we found our typo. Perhaps we should
take a few more steps, before returning to the more general problem.

** Dan**: What does that mean, concretely?

** Carol**: Not bad for one step, perhaps, but our orbit has a radius of unity,
which means a circumference of . With a velocity
of unity, it will take steps to go around the circle,
at the rate we are going. And if every step introduces `only' an
error of , and if the errors built up linearly, we
wind up with a total error of .
That is already a 3% error, even during the first revolution! And after
just a few dozen revolutions, if not earlier, the results will be meaningless.

** Dan**: Good point. Of course, we don't know whether the errors build up
linearly, but for lack of a better idea, that would be the first guess.
Perhaps we should take an even smaller time step. What would happen
if we would use ? Let's repeat your analysis.
After one step, we would have

** Erica**: Bravo! You have just proved that the forward Euler scheme is
a first-order scheme! Remember our discussion at the start? For a
first-order scheme, the errors scale like the first power of the time
step. You just showed that taking a time step that is ten times smaller
leads to a ten times smaller error after completing one revolution.

** Dan**: Great! I thought that numerical analysis was a lot harder.

** Carol**: Believe me, it * is* a lot harder for any numerical integration
scheme that is more complex than first-order. You'll see!

** Dan**: I can wait. For now I'm happy to work with a scheme which I can
completely understand.

** Carol**: So were are we. To sum up: we have verified that our
simple one-step code euler_circular_step.rb does exactly what
we want it to do. And we have validated that what we want it to do is
reasonable: for smaller and smaller time steps the orbit should stay
closer and closer to the true circular orbit.

** Dan**: That's the good news. But at the same time, that's also the bad
news! When we tried to integrate an arbitrary elliptical orbit, we got
a nonsense picture. How come?

** Erica**: We'll have to work our way up to the more complicated situation.
Let us stick to the circular orbit for now. We have a basis to start from:
the first step was correct, that we know for sure. Let's do a few more steps.

** Dan**: Let's try a thousand steps again.

** Carol**: Better to do ten steps. Each time we tried to jump forward too
quickly we've run into problems!

** Erica**: How about a hundred steps, as a compromise? Let's go back
to the very first code we wrote, but now for a circular orbit, and for
an integration of one time unit, which will give us a hundred steps.

** Carol**: Let's keep the old code, for comparison. Here is the new one.
I will call it euler_circular_100_steps.rb:

include Math x = 1 y = 0 z = 0 vx = 0 vy = 1 vz = 0 dt = 0.01 print(x, " ", y, " ", z, " ") print(vx, " ", vy, " ", vz, "\n") 100.times{ r2 = x*x + y*y + z*z r3 = r2 * sqrt(r2) ax = - x / r3 ay = - y / r3 az = - z / r3 x += vx*dt y += vy*dt z += vz*dt vx += ax*dt vy += ay*dt vz += az*dt print(x, " ", y, " ", z, " ") print(vx, " ", vy, " ", vz, "\n") }

** Dan**: If you keep making the names longer and longer, they won't fit on
a single line anymore!

** Carol**: You do what you want and I do what I want; I just happen to sit
behind the keyboard.

** Erica**: Peace, peace! Let's not fight about names; we can later make
copies with shorter names as much as we like.

|gravity> ruby euler_circular_100_steps.rb > euler_circular_100_steps.out

** Dan**: Told you so!

** Carol**: No, you didn't! You didn't give any reason for returning to the
old value.

** Erica**: Hey guys, don't get cranky. Let's just go back to
our original choice of 1,000 steps.

** Carol**: In that case, let's make yet another new file . . . Dan,
close your eyes, I'm adding one more
character to the file name . . . euler_circular_1000_steps.rb:

include Math x = 1 y = 0 z = 0 vx = 0 vy = 1 vz = 0 dt = 0.01 print(x, " ", y, " ", z, " ") print(vx, " ", vy, " ", vz, "\n") 1000.times{ r2 = x*x + y*y + z*z r3 = r2 * sqrt(r2) ax = - x / r3 ay = - y / r3 az = - z / r3 x += vx*dt y += vy*dt z += vz*dt vx += ax*dt vy += ay*dt vz += az*dt print(x, " ", y, " ", z, " ") print(vx, " ", vy, " ", vz, "\n") }

And here is our new result, which I'm calling figure 20:

|gravity> ruby euler_circular_1000_steps.rb > euler_circular_1000_steps.outMuch better!

** Erica**: Indeed, we're really make progress.

** Dan**: We've come around full circle -- almost! At least the particles
are * trying* to orbit each other in a circle, it seems. They're just
making many small errors, that are piling up.

** Erica**: Now that we're getting somewhat believable results, I would like
to make some printouts of our best pictures. Carol, how do you get gnuplot
to make some prints?

** Carol**: That's easy, once you know how to do it, but it is rather
non-intuitive. The easiest way to find out how to do this is to go
into gnuplot and then to type help, and to work your way down the
information about options. To give you a hint, `set terminal`
and `set output`. Of course, if you use gnuplot for the first
time, you would have no way of guessing that * those* are the keywords
you have to ask help for.

** Erica**: That's a general problem with software. Having a help facility
is a good start, but I often find that I need a meta-help facility to
find out how to find out how to ask the right questions to the help
facility. In any case, I'm happy to explore gnuplot more, some day,
but just for now, why don't you make a printout of the last figure we
have just produced.

** Carol**: Okay, here is how it goes. I'm using abbreviations such as
` post` for ` postscript` and ` q` for ` quit`:

|gravity> gnuplot gnuplot> plot "euler_circular_1000_steps.out" gnuplot> set term post eps Terminal type set to 'postscript' Options are 'eps noenhanced monochrome dashed defaultplex "Helvetica" 14' gnuplot> set output "euler_circular_1000_steps.ps" gnuplot> replot gnuplot> qLet's print out this plot:

|gravity> lpr euler_circular_1000_steps.ps

** Carol**: Welcome to the wonderful world of gnuplot. This is a strange
quirk which is one of the things I really don't like about it. But as
long as we use gnuplot, we have to live with it.

Previous | ToC | Up | Next |