10. Multi-Step Methods
Alice: Hi, Bob! You've been up to something again, I can tell.
Bob: I got hooked on exploring multi-step methods.
Alice: Don't tell me that you went to 6th order!
Bob: No, in the sense that I did not just go to 6th order.
I actually went to 6th and 8th order.
Alice: No! That can't be, with your dislike for combinatorics. You're
pulling my leg.
Bob: No, I did not do the combinatorics myself. Instead, I used a
symbolic manipulation package. When I saw you struggling with all those
substitutions when we tracked the fourth-order nature of the
Runge-Kutta integrator, I was already wondering whether we shouldn't
use a symbolic package. At that point I didn't bring it
up, since we were in the middle of trying to figure out how things
worked in the first place, and introducing yet another tool might have
been confusing.
And also when we generalized my two-step method to a four-step method,
I could see the merit of tracking things done with pen and paper.
Certainly when we're going to show this to our students, we should
demand that they are at least able to do it with pen and paper. But
now that the procedure has become crystal clear, we are just faced
with solving a large number of linear equations with the same number
of unknowns, and that is really not something I want to do by hand;
not would I learn much from doing so.
Incidentally, I also had a look at the literature, and I found that
people general solve this type of problem using what they call divided
differentces, really a form of Newton interpolation. In that case
they do not have to solve the system of equations that we set out to
do. Even so, their approach is not that simple either when going to
higher order. Given the ease of using a symbolic manipulation program,
I think that our approach is just as good as the standard one.
Here is what I found, with the results for the coefficients directly
plugged back into the code:
Alice: Again not so bad, and each time derivative still fits on one line.
I see a jerk, a snap, a crackle, and a pop. But what is that next line,
this x variable? It must be the derivative of the pop, given that you
multiply it with yet one higher power of the time step. What does it stand
for?
Bob: To be honest, I had no idea what to call it, so I just called it
x. The advantage was that I could then call the next two derivatives
y and z, for the eighth-order scheme.
Alice: Iron logic! And you used Yoshida's sixth-order integrator, to
rev up the engine, five times.
Bob: Yes, I am glad we had both yo6 and yo8 lying
around. They sure came in handy!
Alice: And you'd like to show me the eighth-order version as well,
don't you?
Bob: I can't wait! Here it is, in file ms8body.rb.
As you can see, this time I abbreviated the variables yet a little
more: instead of @ap3, for example, I wrote @a3,
since by now the meaning is clear, I hope.
Alice: Let's see whether your symbolic manipulator has done its job
correctly . . .
Bob: . . . and whether I have done my job correctly, typing in the
results. We'll start with our, by now almost canonical, set of values:
Alice: That is strange indeed. Independent of the order, the results
should at least get worse, certainly for such a short part of a Kepler
orbit. It's not like we just encountered a singularity or something.
Bob: It feels as if I have just violated the second law of thermodynamics
or something. How can that be?
Alice: As you say in cases like this, let's wiggle the wires, just to
see what happens. How about taking an even larger time step, by a factor
of two?
Bob: May as well. Here is the result:
Bob: If we would have just done the last two runs, we would have been
very happy and content. How difficult it is, to test a piece of software!
We know that there is something seriously wrong, because the first run
of the last three gave a result that was way off, and yet the last two
runs with much larger time steps seem to behave like a charm.
Alice: Like a charm indeed. Now there must be something different,
between the first run, and the last two.
Bob: But what can there be different? Same code, some procedure, same
Bob, same Alice.
Alice: Only larger time steps. And since you kept the total
integration time constant, fewer steps in total. In fact, the last run
had just one step, and the run before that only two steps. But how can
that make a diff . . . . Ah! You talked about revving up the engine . . .
Bob: . . . so if I am only doing one or two time steps, as I've
done in the last two runs . . .
Alice: . . . you're not testing Bob, you're testing Yoshida!
Bob: And clearly, Yoshida has done a very good job indeed. And in the
next-to-last run, where I took ten time steps, I just started to run the
eight-step method, but only in the last three steps. Okay! What a relief.
No problems with the law of thermodynamics, after all!
10.1. An Six-Step Method
10.2. An Eight-Step Method
10.3. Very Strange
|gravity> ruby integrator_driver4g.rb < euler.in
dt = 0.01
dt_dia = 0.1
dt_out = 0.1
dt_end = 0.1
method = ms6
at time t = 0, after 0 steps :
E_kin = 0.125 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 0
(E_tot - E_init) / E_init =-0
at time t = 0.1, after 10 steps :
E_kin = 0.129 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 3.38e-13
(E_tot - E_init) / E_init =-3.86e-13
1.0000000000000000e+00
9.9499478008960474e-01 4.9916426216165405e-02
-1.0020902859905861e-01 4.9748796006154566e-01
Hmm, almost too good. Let's be bold and try a much longer time step,
to get some more playing room, away from machine accuracy.
|gravity> ruby integrator_driver4h.rb < euler.in
dt = 0.05
dt_dia = 0.1
dt_out = 0.1
dt_end = 0.1
method = ms6
at time t = 0, after 0 steps :
E_kin = 0.125 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 0
(E_tot - E_init) / E_init =-0
at time t = 0.1, after 2 steps :
E_kin = 0.129 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 1.53e-13
(E_tot - E_init) / E_init =-1.75e-13
1.0000000000000000e+00
9.9499478009234266e-01 4.9916426209296906e-02
-1.0020902857520513e-01 4.9748796006113827e-01
Hey, that is very strange! Now the errors are getting better,
for a five time longer time step? That can't be!
10.4. Wiggling the Wires
|gravity> ruby integrator_driver4i.rb < euler.in
dt = 0.1
dt_dia = 0.1
dt_out = 0.1
dt_end = 0.1
method = ms6
at time t = 0, after 0 steps :
E_kin = 0.125 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 0
(E_tot - E_init) / E_init =-0
at time t = 0.1, after 1 steps :
E_kin = 0.129 , E_pot = -1 , E_tot = -0.875
E_tot - E_init = 9.12e-12
(E_tot - E_init) / E_init =-1.04e-11
1.0000000000000000e+00
9.9499478026806454e-01 4.9916425775239165e-02
-1.0020902692758932e-01 4.9748796009965129e-01
Alice: And now the error finally gets worse, and by a factor of about 60,
just what you would expect for a 6th-order scheme.
.