RK4 Bouncing a Ball

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?

• Note to moderators/close voters: I did notice the other 'ball bouncing' questions; however all of them are Euler integration. Nov 17, 2011 at 20:38
• Just a quick lookover, it seems like you would bounce 4 times per RK4 call. Nov 17, 2011 at 22:04
• @Jimmy I just tried that :) - it halted the ball (at 167.7*). Nov 17, 2011 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) Nov 17, 2011 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. Nov 18, 2011 at 3:18

• 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? Nov 18, 2011 at 10:58