Code: Select all

```
--Runge-Kutta 4 Lua implementation
--2014 Kyle Emmerich
--
--Reference: http://gafferongames.com/game-physics/integration-basics/
--Special thanks to Glenn Fiedler
local RK4 = {}
--Simple structures for ease of use
RK4.Derivative2D = function(dx, dy, dvx, dvy)
return {dx = dx, dy = dy, dvx = dvx, dvy = dvy}
end
RK4.Derivative1D = function(dx, dv)
return {dx = dx, dv = dv}
end
RK4.State2D = function(x, y, vx, vy)
return {x = x, y = y, vx = vx, vy = vy}
end
RK4.State1D = function(x, v)
return {x = x, v = v}
end
--Evaluators
RK4.Evaluate2D = function(initial, t, dt, deriv, ax, ay)
initial.x = initial.x + deriv.dx * dt
initial.y = initial.y + deriv.dy * dt
initial.vx = initial.vx + deriv.dvx * dt
initial.vy = initial.vy + deriv.dvy * dt
return RK4.Derivative2D(initial.vx, initial.vy, ax, ay)
end
RK4.Evaluate1D = function(initial, t, dt, deriv, a)
initial.x = initial.x + deriv.dx * dt
return RK4.Derivative1D(initial.v, a)
end
--Integrators
RK4.Integrate2D = function(initial, t, dt, ax, ay)
local a = RK4.Evaluate2D(initial, t, 0, RK4.Derivative2D(0, 0, 0, 0), ax, ay)
local b = RK4.Evaluate2D(initial, t, dt * 0.5, a, ax, ay)
local c = RK4.Evaluate2D(initial, t, dt * 0.5, b, ax, ay)
local d = RK4.Evaluate2D(initial, t, dt, c, ax, ay)
local dxdt = 1.0 / 6.0 * (a.dx + 2.0 * (b.dx + c.dx) + d.dx)
local dydt = 1.0 / 6.0 * (a.dy + 2.0 * (b.dy + c.dy) + d.dy)
local dvxdt = 1.0 / 6.0 * (a.dvx + 2.0 * (b.dvx + c.dvx) + d.dvx)
local dvydt = 1.0 / 6.0 * (a.dvy + 2.0 * (b.dvy + c.dvy) + d.dvy)
return RK4.State2D(
initial.x + dxdt * dt,
initial.y + dydt * dt,
initial.vx + dvxdt * dt,
initial.vy + dvydt * dt
)
end
RK4.Integrate1D = function(initial, t, dt, acc)
local a = RK4.Evaluate1D(initial, t, 0, RK4.Derivative1D(0, 0), acc)
local b = RK4.Evaluate1D(initial, t, dt * 0.5, a, acc)
local c = RK4.Evaluate1D(initial, t, dt * 0.5, b, acc)
local d = RK4.Evaluate1D(initial, t, dt, c, acc)
local dxdt = 1.0 / 6.0 * (a.dx + 2.0 * (b.dx + c.dx) + d.dx)
local dvdt = 1.0 / 6.0 * (a.dv + 2.0 * (b.dv + c.dv) + d.dv)
return RK4.State1D(
initial.x + dxdt * dt,
initial.v + dvdt * dt
)
end
return RK4
```