Take the 2-minute tour ×
Game Development Stack Exchange is a question and answer site for professional and independent game developers. It's 100% free, no registration required.

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?

share|improve this question
    
Note to moderators/close voters: I did notice the other 'ball bouncing' questions; however all of them are Euler integration. –  Jonathan Dickinson Nov 17 '11 at 20:38
    
Just a quick lookover, it seems like you would bounce 4 times per RK4 call. –  Jimmy Nov 17 '11 at 22:04
    
@Jimmy I just tried that :) - it halted the ball (at 167.7*). –  Jonathan Dickinson Nov 17 '11 at 22:09
    
I think the acceleration shouldn't be based on t, but perhaps modeled as some kind of energy-conserving spring force K*(state.X - 100) –  Jimmy Nov 17 '11 at 23:27
    
If you're doing a bouncing ball in gravity then acceleration will be a constant -9.8m/s^2 to represent the downwards pull. The velocity of the ball changes based on acceleration and time. The position of the ball changes based on velocity and time. If/When the ball hits the floor (let's say 0.0f) the ball's velocity will reflect back, so velocity = -velocity when that happens. The article you linked has a spring with a damper, a more complicated set of accelerations than your bouncing ball in gravity. –  Patrick Hughes Nov 18 '11 at 3:18

1 Answer 1

up vote 1 down vote accepted

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.

share|improve this answer
    
I'll give that a shot tonight. I guess at that point I need to zero the accumulated time and set the initial state (Evaluate(eulerIntegrator, ref state, t, 0.0f, ***ref a***, ref a)) to the state when the bounce occurred? –  Jonathan Dickinson Nov 18 '11 at 10:58

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.