./ch07.html

ToC

7.1. One Integration Step: Verification

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

Dan: . . . and getting just one short new line of output, after the initial conditions are echoed. Nice and simple!

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

and

so the new position must be:

and since we are using

, we expect:

Carol: And this is indeed what we see in the first half of the second output line.

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:

In our case this gives

And this in turn means that

Carol: And this is not what is printed in the last half of the last output line, in our one-step program.

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

Dan: Wonderful! Now both the position and the velocity components are correct, after the first step. We are winning!

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

7.2. A Different Surprise

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.out

And here is the plot, in fig 18:

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

Figure 18: Second attempt at integrating the two-body problem: a different failure.

Carol: Well, Dan, what do you say?

Dan: No comment.

7.3. One Integration Step: Validation

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.

Carol: I agree. Let us sum up what we've done with our one-step code. We have verified that our program does what the algorithm intended, and that is certainly nice! But it is only half of the work. We now have to check whether our particular algorithm does indeed give a reasonable result, which corresponds to the behavior of gravitating particles in the real world. This is called validation. In the computer science literature these two checks are often called V&V, for Verification and Validation.

In other words, so far in our one-step program we have passed the verification test. The computer code does exactly what we wanted it to do, at least for that one step. But now we have to do a validation test.

Dan: What does that mean, concretely?

Carol: For example, we can ask whether the first step keeps the two particles on a circular orbit. We can answer that question with pure thought. After the first step, the new separation is:

Dan: Instead of the correct value of

, we are half a one hundredth of one percent off. Not bad, I would say.

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

We now need roughly

steps to go around the circle, If the errors build up linearly, the radial separation will grow to something like

. Aha! Only a 0.3% error, instead of 3%.

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.

7.4. More Integration Steps

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

Figure 19: Third attempt at integrating the two-body problem: part of a circle?

Carol: Figure 19 may or may not be part of a good circle; hard to see when you only have a small slice.

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.

7.5. Even More Integration 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.out

Figure 20: Fourth attempt at integrating the two-body problem: looking much better.

Much 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.

7.6. Printing Plots

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> q

Let's print out this plot:

  |gravity> lpr euler_circular_1000_steps.ps

Erica: Great, same figure. But hey, wait a minute, the symbol used for the points in the printed figure is very different from the symbol that appeared on the screen!