Previous ToC Up Next

17. Modified Euler in Vector Form

17.1. An Easy Translation

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:

In our first code this became: 

   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:

Dan: It was a lot of work, but now that we got the vector class, I must admit that the code looks a lot more readable. So I guess this will make life a lot easier for us. But before we move on, does it give the correct answers?

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.0
and 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

17.2. Variable Names

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

17.3. Consistency

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
Dan: Sure, that is more consistent, but you've just made the code two lines longer! In fact, five lines longer, if you count the three blank lines. Why did you add blank lines?

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!

17.4. A Method Returning Multiple Values

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
Dan: I'm glad it works, but how does it work? It seems that your new method step_pos_vel returns an array.

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!

17.5. Simplification

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:

Figure 35: Two-step implementation of the modified Euler algorithm

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:

Figure 36: One-step implementation of the modified Euler algorithm

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
Carol: And so it does, almost. Since some of the additions and multiplications and such are done in a different order, round-off errors may prevent us from getting the exact same results, but this is certainly close enough. And yes, I admit, this code is quite a bit simpler.

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