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?
