RK4 Bouncing a Ball

Question

I am trying to wrap my head around RK4. I decided to do the most basic 'ball with gravity that bounces' simulation. I have implemented the following integrator given Glenn Fiedler's tutorial:

/// <summary>
/// Represents physics state.
/// </summary>
public struct State
{
    // Also used internally as derivative.

    // S: Position
    // D: Velocity.
    /// <summary>
    /// Gets or sets the Position.
    /// </summary>
    public Vector2 X;

    // S: Position
    // D: Acceleration.
    /// <summary>
    /// Gets or sets the Velocity.
    /// </summary>
    public Vector2 V;
}

/// <summary>
/// Calculates the force given the specified state.
/// </summary>
/// <param name="state">The state.</param>
/// <param name="t">The time.</param>
/// <param name="acceleration">The value that should be updated with the acceleration.</param>
public delegate void EulerIntegrator(ref State state, float t, ref Vector2 acceleration);

/// <summary>
/// Represents the RK4 Integrator.
/// </summary>
public static class RK4
{
    private const float OneSixth = 1.0f / 6.0f;

    private static void Evaluate(EulerIntegrator integrator, ref State initial, float t, float dt, ref State derivative, ref State output)
    {
        var state = new State();

        // These are a premature optimization. I like premature optimization.
        // So let's not concentrate on that.
        state.X.X = initial.X.X + derivative.X.X * dt;
        state.X.Y = initial.X.Y + derivative.X.Y * dt;
        state.V.X = initial.V.X + derivative.V.X * dt;
        state.V.Y = initial.V.Y + derivative.V.Y * dt;

        output = new State();
        output.X.X = state.V.X;
        output.X.Y = state.V.Y;
        integrator(ref state, t + dt, ref output.V);
    }

    /// <summary>
    /// Performs RK4 integration over the specified state.
    /// </summary>
    /// <param name="eulerIntegrator">The euler integrator.</param>
    /// <param name="state">The state.</param>
    /// <param name="t">The t.</param>
    /// <param name="dt">The dt.</param>
    public static void Integrate(EulerIntegrator eulerIntegrator, ref State state, float t, float dt)
    {
        var a = new State();
        var b = new State();
        var c = new State();
        var d = new State();

        Evaluate(eulerIntegrator, ref state, t, 0.0f, ref a, ref a);
        Evaluate(eulerIntegrator, ref state, t + dt * 0.5f, dt * 0.5f, ref a, ref b);
        Evaluate(eulerIntegrator, ref state, t + dt * 0.5f, dt * 0.5f, ref b, ref c);
        Evaluate(eulerIntegrator, ref state, t + dt, dt, ref c, ref d);

        a.X.X = OneSixth * (a.X.X + 2.0f * (b.X.X + c.X.X) + d.X.X);
        a.X.Y = OneSixth * (a.X.Y + 2.0f * (b.X.Y + c.X.Y) + d.X.Y);
        a.V.X = OneSixth * (a.V.X + 2.0f * (b.V.X + c.V.X) + d.V.X);
        a.V.Y = OneSixth * (a.V.Y + 2.0f * (b.V.Y + c.V.Y) + d.V.Y);

        state.X.X = state.X.X + a.X.X * dt;
        state.X.Y = state.X.Y + a.X.Y * dt;
        state.V.X = state.V.X + a.V.X * dt;
        state.V.Y = state.V.Y + a.V.Y * dt;
    }
}

After reading over the tutorial I noticed a few things that just seemed 'out' to me. Notably how the entire simulation revolves around t at 0 and state at 0 - considering that we are working out a curve over the duration it seems logical that RK4 wouldn't be able to handle this simple scenario. Never-the-less I forged on and wrote a very simple Euler integrator:

static void Integrator(ref State state, float t, ref Vector2 acceleration)
{
    if (state.X.Y > 100 && state.V.Y > 0)
    {
        // Bounce vertically.
        acceleration.Y = -state.V.Y * t;
    }
    else
    {
        acceleration.Y = 9.8f;
    }
}

I then ran the code against a simple fixed-time step loop and this is what I got:

0.05 0.20 0.44 0.78 1.23 1.76 ... 74.53 78.40 82.37 86.44 90.60 94.86 99.23 103.05 105.45 106.94 107.86 108.42 108.76 108.96 109.08 109.15 109.19 109.21 109.23 109.23 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 109.24 ...

As I said, I was expecting it to break - however I am unsure of how to fix it. I am currently looking into keeping the previous state and time, and working from that - although at the same time I assume that will defeat the purpose of RK4.

How would I get this simulation to print the expected results?

Answer 1

Try this: run your Runge-Kutta integration to compute the motion, and then when it's done check for the bounce condition and update your state. It worked for me. And yes, this does mean that you still have to patch up the inter-penetration like you would with Euler integration.

I can't claim this authoritatively, but I think step changes don't go well with the RK4 state evaluator function.