Accuracy vs. timestep for different techniques

Let's use the "satellite motion around the Earth" project to test the accuracy of several different techniques for computing the motion. We'll include

• E1 Euler's method to first order. We compute the new position like so:
  

• E1a The "symplectic" Euler's method to first order. In this case, we use the new velocity to compute the new position. A simple change, but it increases the accuracy quite a bit.
  

• E2 Euler's method to second order.
  

• RK4 Fourth-order Runge-Kutta method, in which we first predict a set of future positions at both one-half and one full timestep in the future, and use those predictions to compute several velocity-like quantites k₁ to k₄. We can then compute the new position with a weighted average:
  

To test these different methods, we'll follow the motion of a satellite around the Earth. We assume a nicely spherical Earth, ignore air resistance and all other complications, and begin with initial conditions:

• satellite at an altitude of H = 250 km above the surface
• satellite moving with velocity v₀ = 7762 m/s, which ought to produce a perfectly circular orbit

We'll follow the motion of the satellite for 10 days = 864,000 seconds, which corresponds to roughly 160 orbits. Before we start, we'll compute the satellite's initial total energy E₀

At the end of the simulation, we compute the final total energy E_f. In a perfect simulation, the initial and final energies would be exactly the same. In the actual simulation, they aren't, due to the accumulated errors in each step of the motion. We use the fractional change in total energy f

as a metric to describe the quality of the simulation. The smaller f is, the better.

Now, to investigate the properties of each technique, we'll run the simulation many times, using a fixed timestep in each trial. For each technique, we will run trials with timesteps between 0.0001 and 1000 seconds, doubling the timestep each time. So, for example,

• first trial has timestep 0.0001 seconds
• second trial has timestep 0.0002 seconds
• third trial has timestep 0.0004 seconds
...
• 23'rd trial has timestep 419 seconds
• 24'rd trial has timestep 838 seconds

  Q:  How many steps does the first trial take?  
      Into how many pieces is each orbit broken?

  Q:  How many steps does the final trial take?  
      Into how many pieces is each orbit broken?

The answers

After running all these simulations (it took about 28 hours on an Intel Pentium 4 3GHz CPU), one can compare the fractional error as a function of timestep.

The results have a number of interesting features.

• When the fractional change f is close to 1, then the simulation is basically worthless. Every method reaches this point when the timestep grows too large, but the higher-order methods have much larger critical timesteps.
• Two of the methods, E2 and RK4, appear to grow WORSE with smaller timesteps, for step sizes below some value. Can you explain why?
• Every method does show a linear pattern over some wide range of timesteps. Since we are plotting logarithmic quantities here, that means that within those ranges, there's a relationship like this:
  

  But we can also express the relationship in a more convenient way:

  

  Q:  What are the slopes for each technique
      over the linear portions of their results?




  Q:  For each technique, complete the following
      statement:

        If one cuts the timestep in half,
          the technique becomes ______
          times as accurate.

When the timesteps are very large, none of the methods works very well. Even though the change in fractional energy may be small, a simulation may not reflect reality. Conservation of energy does not guarantee correct results! Look at the first "orbit" for each technique if we use a timestep of 500 seconds:

Here's a more detailed view, showing the outline of the Earth.

The symplectic Euler's method had a value of f = 0.1 -- yet it sends the satellite off into a spiral trajectory. The second-order Euler's method claims to have f = 0.04, but it is obviously even worse. The Runge-Kutta method yields f = 0.01, which is almost the same as that for second-order Euler ... but its trajectory is ALMOST a perfect circle.

Don't fall into the trap of trusting the numbers. Always examine results with the Mark I eyeball -- it's an excellent detector of bogus results.