3.4 Extending a Code

Carol:

In addition, rather than changing the program each time we change step size, I'll make it an option to be provided by the user. How about this version:

Bob:

And even with a polite request to remind the user what is needed. Very nice. Let's compile, run, and plot!

Carol:

Okay. Let's direct to output to a file forward2a.out, and let's rename the old output file forward1.out, to get a consistent numbering system, where the number in the source file corresponds to the same number in the output file; the `a' after `2' is added here since different input options will produce different output files.

Bob:

How about making it even clearer and add the value of the time step that we input to the name?

Carol:

Yes, even better. Let's start with $dt = 0.01$, the same value we used by default in forward_euler1.C. It never hurts to check whether we indeed get the same result. Here goes:

|gravity> g++ -o forward_euler2 forward_euler2.C
|gravity> mv forward.out forward1.out
|gravity> forward_euler2 > forward2_0.01.out

Bob:

So much for the polite reminder. What is happening, is the computer so slow? Or is it just hanging?

Alice:

Ah, I see. By redirecting the output to forward2_0.01.out, we are also redirecting our reminder to that file, spoiling the format of our nbody snapshot output, and losing our reminder.

Carol:

Yes, of course, you are right. Well, how about putting the reminder on the error stream cerr, rather than the output stream cout? While reminding you of something is not an error, it does provide a handy way to disentangle the snapshot output, meant to be read further by a computer, and the reminder output, aimed at human consumption.

Bob:

Good idea. But while we do that, let's make another modification as well. I just realized that our halting criterion is completely wrong. If we specify dt = 0.0001, the system would bail out after only 1000 steps, as before, thus advancing the time from zero to 0.1, barely enough to see the particles move, and far too little to complete even one orbit.

Alice:

Right you are. Let's put in a time counter, so that we can keep integrating till time 10, for good measure.

Carol:

Fine! And I'll rename the previous version forward_euler2a.C, while keeping the name forward_euler2.C for the corrected version. When you are developing code, it is always a good idea to keep older versions lying around, even if they are wrong. After all, subsequent versions may turn out to be even more wrong, in which case you might well want to go back to the less wrong one! Here we are:

Carol:

And now let's try again:

Bob:

That's curious. Almost the same results as before, but not quite. Remember, earlier we got:

Bob:

Ah, I see: same results, after all, but one more time step. Our next to last line is identical to our previous last line.

Carol:

That would imply that the new file should be one line longer. Let's check by doing a word count (actually a line-word-character count):

Alice:

Indeed: one line more and six words more, for the extra position and velocity components. I guess our switch from integer book keeping with the number of time steps, to floating-point book keeping by using the time introduced some slight round-off which caused the program to do one more step.

Bob:

Sounds plausible, but let's check it out. Seeing is believing! Shall we put in an extra output statement on cerr? This time we're after errors after all, namely round-off errors. Let's add the following line at the end, after the two lines starting with cout:

	cerr << t << endl;

Carol:

Done. I'll call the resulting code forward_euler2b.C, following our tradition to keep all version around for a while. Let's run the code again, but I won't bother about the snapshot output, since we already know what that is. Redirecting the cout output to /dev/null is a time-honored Unix way to get rid of it. The `null' device acts like a black hole, getting rid of whatever falls into it, without a trace. We are then left only with the cerr output on the screen.

Bob:

Hah! So much for your clever theory, Alice. It stopped at t = 10, right on the mark. So round-off cannot have been the reason.

Alice:

Not so quick in claiming victory, Bob. Look at the code. In the for loop the halting criterion is t < 10. If there would have been no round-off, the code should have stopped at t = 9.99, not at t = 10.

Bob:

That makes sense. And yet there is no round-off. It seems like we are both right, but that is impossible. What gives?

Carol:

In case of doubt, find it out. I have a sneaky suspicion that there is another round-off that is fooling us: the limited accuracy of the output. The quickest way to find out is to modify the time output line. Instead of simply printing the time, let's print the difference of the real time t and the non-round-off value t = 10. So our last error output line becomes:

	cerr << t - 10.0 << endl;

Bob:

Okay, Alice was right after all. Close but no cigar. After 1,000 steps, the time is updated to a hair's breadth smaller than 10, and the code takes one extra step. Clever detective work, Carol!

Carol:

I'm glad my course in numerical methods was good for something!

Bob:

By the way, returning to the commands you gave earlier, just before we spotted the extra output line, what was that all about, that you had to remove files by hand?

Carol:

The operating system refused, fortunately, to override the existing output file that we created with our earlier version of forward_euler2.C. This is a good thing, since it is all to easy to override (much) earlier results that are perhaps hard to recreate. Similarly, when I gave the explicit command rm to remove the file, the system checked to make sure I knew what I was doing, and I answered y for yes.

Bob:

That seems like overkill to me. You asked for a specific file to be deleted, and it checked to see whether you had learned to read and write??

Carol:

For a single file this may be overkill, but often we give a command such as rm tmp* if we want to remove all temporary files in a directory. If you (or your cat) were to hit the space bar by mistake in the middle, you could wind up with rm tmp * ...

Bob:

...I see, and thus deleting all files in that directory. No, that would not be good, I agree.

Carol:

By the way, you can choose this safety option by typing rm -i with option i for inquiry, prompting the system to question you for each file to be deleted. Or simpler, you can add this to your shell definitions, for example to your .cshrc file if you use a C shell.

Bob:

And what about !$?

Carol:

Oh, that is shorthand for retyping the last word in the previous command. I could have used that two lines earlier too, by typing rm !$. It speeds things up, just like the command !! which repeats the whole previous command. Also the second time I ran the program I could have simply typed !f which would have expanded to forward_euler2 > forward2_0.01.out. Or I could have type !-2 to repeat the command that was issued two steps before. Unix has many ways to say a lot with two or three key strokes! Here is what I could have typed: