Previous | ToC | Up | Next |

** Dan**: Now that we have a vector version of forward Euler, it's time
to clean up our modified Euler code as well.

** Carol**: That will be an easy translation. I will start by copying
the old code from euler_modified_array.rb into a new file,
euler_modified_vector_try1.rb. All we have to do is to translate
the code from array notation to vector notation. Same as what we did
with euler_vector.rb. Here it is:

require "vector.rb" include Math def print_pos_vel(r,v) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} print "\n" end r = [1, 0, 0].to_v v = [0, 0.5, 0].to_v dt = 0.01 print_pos_vel(r,v) 1000.times{ r2 = r*r r3 = r2 * sqrt(r2) a = -r/r3 r1 = r + v*dt v1 = v + a*dt r12 = r1*r1 r13 = r12 * sqrt(r12) a1 = -r1/r13 r2 = r1 + v1*dt v2 = v1 + a1*dt r = 0.5 * ( r + r2 ) v = 0.5 * ( v + v2 ) print_pos_vel(r,v) }

** Erica**: What a relief! The lines are shorter,
there are fewer lines, but what is most important: the lines are easy
to understand, with a direct correspondence between code and math.

Let's trace our history, in this regard. We started off writing with pen and paper:

x1 = x + vx*dt y1 = y + vy*dt z1 = z + vz*dt

Then in our array code it became

r1 = [] r.each_index{|k| r1[k] = r[k] + v[k]*dt}

and finally, in our vector code, we wrote:

r1 = r + v*dt

which is very close indeed to what we started out with:

** Carol**: Here's the old result, from the array code:

|gravity> ruby euler_modified_array.rb | tail -1 0.400020239524913 0.343214474344616 0.0 -1.48390077762002 -0.0155803976141248 0.0and here's what our new vector code gives:

|gravity> ruby euler_modified_vector_try1.rb | tail -1 0.400020239524913 0.343214474344616 0.0 -1.48390077762002 -0.0155803976141248 0.0

** Dan**: Good! I'm happy.

** Erica**: But I'm not, at least not completely. Look, in the code
we are using the variable `r2` in two very different ways.
Early on, we use it to hold the value of ,
the square of the original variable , defined as
the inner product
. But later, toward the end
of the loop, we use the same variable to hold value of
, the predicted value of .

I guess Ruby doesn't mind that we assign completely different values, even with different types, first a scalar, then a vector. But I sure do mind! And someone else reading our code from scratch is likely to be confused.

** Carol**: You have a point there. Okay, how about calling the initial
position instead of ? That is more
consistent anyway. We can then use the variable name `r0`
instead of `r` for the initial vector, and the scalar value of
its square will then become `r02`. So there will be no possible
confusion anymore! Here is the new listing, in file
euler_modified_vector_try2.rb:

require "vector.rb" include Math def print_pos_vel(r,v) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} print "\n" end r = [1, 0, 0].to_v v = [0, 0.5, 0].to_v dt = 0.01 print_pos_vel(r,v) 1000.times{ r02 = r*r r03 = r02 * sqrt(r02) a = -r/r03 r1 = r + v*dt v1 = v + a*dt r12 = r1*r1 r13 = r12 * sqrt(r12) a1 = -r1/r13 r2 = r1 + v1*dt v2 = v1 + a1*dt r = 0.5 * ( r + r2 ) v = 0.5 * ( v + v2 ) print_pos_vel(r,v) }

and here is the test:

|gravity> ruby euler_modified_vector_try2.rb | tail -1 0.400020239524913 0.343214474344616 0.0 -1.48390077762002 -0.0155803976141248 0.0

** Erica**: Yes, that's better, and you are no longer using the same
variable name for two different things. But you haven't quite done
what you said you would do, namely calling the initial position
instead of . You have only assigned
the square of to , or `r02`
in the code.

** Carol**: That's because I wanted to continue using the original variables
` r` and ` v`, to keep track of the evolving code. The alternative would
have been to call the running variables `r0` and `v0`, but
that would be misleading, as if the particle would come back to the original
position each time.

** Erica**: How about a compromise? We can keep the original variables
` r` and ` v`, but convert them to `r0` and `v0` at the
beginning of the loop. Let me try this, in file
euler_modified_vector_try3.rb:

require "vector.rb" include Math def print_pos_vel(r,v) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} print "\n" end r = [1, 0, 0].to_v v = [0, 0.5, 0].to_v dt = 0.01 print_pos_vel(r,v) 1000.times{ r0 = r v0 = v r02 = r0*r0 r03 = r02 * sqrt(r02) a0 = -r0/r03 r1 = r0 + v0*dt v1 = v0 + a0*dt r12 = r1*r1 r13 = r12 * sqrt(r12) a1 = -r1/r13 r2 = r1 + v1*dt v2 = v1 + a1*dt r = 0.5 * ( r0 + r2 ) v = 0.5 * ( v0 + v2 ) print_pos_vel(r,v) }

and let me test it right away:

|gravity> ruby euler_modified_vector_try2.rb | tail -1 0.400020239524913 0.343214474344616 0.0 -1.48390077762002 -0.0155803976141248 0.0

** Erica**: I'm not really worried about the extra two lines of code.
What's much more important is that in this new notation we can see
clearly that we have a new target of attack for the DRY principle.
Look, the two blocks of code, that I have highlighted by place blank
lines around them, are nearly identical!

** Carol**: Right you are! This calls for a new method. Here, let me try,
in file euler_modified_vector.rb:

require "vector.rb" include Math def print_pos_vel(r,v) r.each{|x| print(x, " ")} v.each{|x| print(x, " ")} print "\n" end def step_pos_vel(r,v,dt) r2 = r*r r3 = r2 * sqrt(r2) a = -r/r3 [r + v*dt, v + a*dt] end r = [1, 0, 0].to_v v = [0, 0.5, 0].to_v dt = 0.01 print_pos_vel(r,v) 1000.times{ r1, v1 = step_pos_vel(r,v,dt) r2, v2 = step_pos_vel(r1,v1,dt) r = 0.5 * ( r + r2 ) v = 0.5 * ( v + v2 ) print_pos_vel(r,v) }

and I'll test it right away:

|gravity> ruby euler_modified_vector.rb | tail -1 0.400020239524913 0.343214474344616 0.0 -1.48390077762002 -0.0155803976141248 0.0

** Carol**: Ah, yes. In Ruby you can return more than one value simultaneously,
and the way to do that is to list them in an array. Then, when you invoke
the method, you can just list the variables to which you want to assign the
array values, in the same order. Very lazy, very simple, and in a way, you
might say, following the principle of least surprise.

** Dan**: Well, yes, once you realize what is happening.

** Erica**: And it is quite elegant, I must say. I begin to like Ruby more and
more! The inner loop now has become almost a mathematical formula, rather
than a piece of computer code. You can read it aloud: step forward, step
again, and then average the values from the beginning and the end, both
for position and velocity, and print the results. Wonderful!

** Erica**: Before we move on, there is something that bothers me, which I
noticed as soon as we translated the modified Euler scheme into vector
notion, in euler_modified_vector_try1.rb. I had not noticed
any problem when we first wrote the much longer component-based version,
but when it became so compact, I realized that we are being clumsy.

** Carol**: How so?

** Erica**: Look, we have effectively implemented figure 35,
right? Let me sketch it here again:

We take two steps forward, in order to compute a single improved step.

However, we started off that whole discussion, way back when, we the simpler figure 36. Let me draw that one again too:

You see, in that original figure, we average two steps, in order to compute an improved step.

** Carol**: But you needed to compute the upper step, before you could compute
the lower step, so really both ways of drawing the picture involve two steps.

** Erica**: True, but to me at least the original figure suggests a somewhat
simpler procedure.

** Dan**: I don't care to quibble about philosophy of aesthetics.
Here is the inner loop of the first vector version that we wrote for
modified Euler:

1000.times{ r2 = r*r r3 = r2 * sqrt(r2) a = -r/r3 r1 = r + v*dt v1 = v + a*dt r12 = r1*r1 r13 = r12 * sqrt(r12) a1 = -r1/r13 r2 = r1 + v1*dt v2 = v1 + a1*dt r = 0.5 * ( r + r2 ) v = 0.5 * ( v + v2 ) print_pos_vel(r,v) }

Where do you think you can simplify things?

** Erica**: Let me copy that file, euler_modified_vector_try1.rb,
to euler_modified_vector_try4.rb, and let me see whether I can
change the loop, to make it look more like my understanding of the
original figure, the second one above:

1000.times{ r2 = r*r r3 = r2 * sqrt(r2) a = -r/r3 r1 = r + v*dt r12 = r1*r1 r13 = r12 * sqrt(r12) a1 = -r1/r13 r += v*dt + 0.5*a*dt*dt v += 0.5*(a + a1)*dt print_pos_vel(r,v) }

** Dan**: That sure looks a lot simpler. Does it give the same answer?

|gravity> ruby euler_modified_vector_try4.rb | tail -1 0.400020239524793 0.34321447434461 0.0 -1.48390077762024 -0.0155803976143043 0.0

** Erica**: Not only that, you can now understand the mathematical
structure better. The increments in position and velocity, in
the last two lines, are just Taylor series expansions, up to terms
that contain the accelerations. In the case of the position, the
acceleration term is second order in ` dt` and so the original value
of the acceleration is good enough. In the case of the velocity,
we need more precision, so we take the average of the forward and
backward Euler values.

** Carol**: On the other hand, the inner loop in our code in
euler_modified_vector.rb was even shorter. Here it is:

1000.times{ r1, v1 = step_pos_vel(r,v,dt) r2, v2 = step_pos_vel(r1,v1,dt) r = 0.5 * ( r + r2 ) v = 0.5 * ( v + v2 ) print_pos_vel(r,v) }

** Dan** Yes, but at the expense of introducing an extra method, making
the whole code longer again.

** Erica**: I guess there is a trade off here. Introducing an extra
method gives the elegance of seeing two steps being taken explicitly,
which shortening the code as I just did brings out the Taylor series
character of the result of averaging two steps.

** Carol**: I agree. Well, we've learned a lot, and I am actually happy
to have several versions. In general, there is never just one right
solution for any question concerning software writing!

Previous | ToC | Up | Next |