Simulation Class: The Euler-Cromer Time Step

Forward Euler Problems

In the previous class, we observed that a Forward Euler time step quickly produced an unusable solution to our simple spring system. This illustration from the Wikipedia entry on Forward Euler illustrates how the method works, and also why its solutions drift away from the correct result.

In this illustration, you can see that at each point, the slope of the function is used to construct a line segment that estimates the position of the next point. The error in each estimate is accumulated with each subsequent estimate, and the solution veers off course. In order to study these systems better, we first need some debug tools.

Adding Debug Tools to our Spring Simulation

The Forward Euler Time Step that we added to the Spring Simulation in the previous class produced an oscillatory motion, but did so in a way that produced an unstable result, with the mass oscillating wildly off screen in a short period of time.

The instability shown by our simple spring simulation is a constant problem for most simulations. It's easy to see the problem in this simple example (which is why we're using it!), but the problem affects large-scale fluid and cloth sims, rigid body sims, and so on. Strangely, despite the general knowledge that Forward Euler is not a good time integration tool, it is still used in many commercial tools, including Maya and Houdini.

In order to analyze our improvements, we first need a measuring stick. To quote my good friend Adam Shand, "You can't improve what you can't measure.". We chose the simple spring because it did have an exact answer. Therefore, we'll add some code to our drawing routine that draws the right answer just above the simulated spring. Here's that code - note that this uses the expression for the explicit position of the spring that we established in the previous class. We'll also add some text to the top of the drawn image to remind us which time integration we're using.

void DrawState()
{
    // Compute end of arm.
    float SpringEndX = PixelsPerMeter * State[StatePositionX];

    // Compute the CORRECT position.
    float sqrtKoverM = sqrt( Stiffness / BobMass );
    float x0 = InitState[StatePositionX];
    float v0 = InitState[StateVelocityX];
    float t = State[StateCurrentTime];
    float CorrectPositionX = ( x0 * cos( sqrtKoverM * t ) ) +
        ( ( v0 / sqrtKoverM ) * sin( sqrtKoverM + t ) );
    
    // Compute draw pos for "correct"
    float CorrectEndX = PixelsPerMeter * State[CorrectPositionX];

    // Draw the spring.
    strokeWeight( 1.0 );
    line( 0.0, 0.0, SpringEndX, 0.0 );
          
    // Draw the spring pivot
    fill( 0.0 );
    ellipse( 0.0, 0.0, 
             PixelsPerMeter * 0.03, 
             PixelsPerMeter * 0.03 );
    
    // Draw the spring bob
    fill( 1.0, 0.0, 0.0 );
    ellipse( SpringEndX, 0.0, 
             PixelsPerMeter * 0.1, 
             PixelsPerMeter * 0.1 );

    // Draw the correct bob in blue
    fill( 0.0, 0.0, 1.0 );
    ellips( CorrectEndX, -PixelsPerMeter * 0.25,
            PixelsPerMeter * 0.1,
            PixelsPerMeter * 0.1 );
}

And here that is:

I also added a "reset" capability to the simple spring code by using Processing's keyRelease function. In response to the key 'r' being released in the view (which has ascii-code 114), we copy the initial state onto the state. (See how handy state vectors are turning out to be!) Because we included the Current Time in our state vector, this is easy:

// Reset function. If the key 'r' is released in the display, 
// copy the initial state to the state.
void keyReleased()
{
    if ( key == 114 )
    {
        for ( int i = 0; i < StateSize; ++i )
        {
            State[i] = InitState[i];
        }
    }  
}

The Euler-Cromer Time Step

When we created the forward finite difference approximation to our ODE equation of motion in the previous class, we had to create two forward difference equations - one for velocity from acceleration, and then one for position from velocity. As you remember, we computed acceleration from the current position, then updated position based on the current velocity, and then finally updated velocity based on the acceleration. What if we updated velocity first, and then updated the position from the current velocity? That would look like this:

\( a( t ) = f( t, x( t ) ) = \frac{-k}{m} x( t ) \\ v( t + \Delta t ) \approx v( t ) + \Delta t a( t ) \\ x( t + \Delta t ) \approx x( t ) + \Delta t v( t + \Delta t ) \)

This seems almost trivially different from our Forward-Euler time step. This method is called many different things: Semi-Implicit Euler Method, Symplectic Euler, Semi-Explicit Euler, Euler-Cromer, and Newton–Størmer–Verlet (NSV) (at least). One of the things I've found in years of writing simulations is that just learning the names of all these different techniques, and learning that the different names often mean exactly the same thing, is a huge part of simplifying the giant mess of understanding.

The code changes to implement Euler Cromer are trivial, we just reorder the time step - literally just swap the velocity update line and the position update line:

// Euler-Cromer Time Step function.
void TimeStep( float i_dt )
{
    // Compute acceleration from current position.
    float A = ( -Stiffness / BobMass ) * State[StatePositionX];

    // Update velocity based on acceleration.
    State[StateVelocityX] += i_dt * A;

    // Update position based on current velocity.
    State[StatePositionX] += i_dt * State[StateVelocityX];

    // Update current time.
    State[StateCurrentTime] += i_dt;
}

And let's see how this ONE LINE (okay, TWO-LINE) change affects the sim:

This is not bad! Our simulation seems to have become stable, which is good news. However, if we watch the simulation for a while, it becomes clear that the simulated positions diverge more and more over time from the correct positions. In order to fix this, we're going to need to explore even higher-order time integration techniques, which is the topic of the next class!

Download

All of the files associated with this class are available via a public github repository, found here: https://github.com/blackencino/SimClass