RK4 Bouncing a Ball

Posted by Jonathan Dickinson on Game Development See other posts from Game Development or by Jonathan Dickinson
Published on 2011-11-17T20:37:31Z Indexed on 2011/11/18 2:04 UTC
Read the original article Hit count: 756

Filed under:
|

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?

© Game Development or respective owner

Related posts about c#

Related posts about physics