Runge-Kutta Integration and Solution for Physical Modelling in an SSBO-based Boid Simulator
Notes by Kate Goss, for panda3d-ssbo
Using Runge-Kutta 4th order helps account for change in velocity over a timestep. It can be extended to offer similar returns for force. Through estimating the rate of change of the independent variable four times over a given change in the dependent variable (e.g. a timestep), it calculates the area under the graph as if the variables change continuously. In the naive case, the speed of an entity changes whenever it is set, so it is displaced by the previous velocity until the next is determined, and so on. In 'real life', the velocity would gradually change, meaning that the position goes up or down faster or slower across the step. RK4 emulates that gradual change from one value to another, giving more realistically physical change in position.
The process can be done by calling a function, the Ordinary Differential Equation (ODE), with the values needed to calculate the new values (inferring a rate of change, not expressed). It then returns a result, which is used to calculate the next result in the chain by suggesting what it would transform to under the initial conditions in half as much time, and this new result is used to repeat this task with a rate likely to be off by an opposing degree to some of the existing error. Using this new rate to calculate the change by the end of the timestep gives us our fourth value, which can be used to create the (1 + 2 + 2 + 1)/6 weighted average.
The pseudocode for the simpler case shows a thorough example that takes a weighted average of both position and velocity to demonstrate another level of additional 'smoothing' that can be achieved using an integrator like RK. It ensures that both the change in acceleration and the 'consequent' change in velocity across the timestep happen gradually, rather than being held constant across the whole timestep. Note, in physical systems, acceleration doesn't 'cause' velocity to change, it is the rate of change of velocity- however, in this computation, we must assign a value to our velocity manually, and therefore the variable representing its rate of change can only act in that capacity by assigning a new value. This is why we modify the value with a weighted average; as if it were reflecting an ongoing mutation in the value, not merely forcing a new state onto it. Extra efficiency is harvested here by running two levels of integration (accel > vel, vel > pos) in one integrator, allowing both estimates to influence each other throughout the timestep, resulting in a more accurate sum for pos at the end.
In order to properly model changes in force over the timestep, a shader must be triggered to run
the ode for each stage of rk4; this will allow the program to access updated positions of
neighbours from the SSBO, appropriately culling calculations with the spatial hashing algorithm.
This methodology works to recognise that the boids would experience continuous forces much like
their velocities change continuously, calculating a smoother average of the effects across the
timestep. Although this runs the calculation 4x per frame, this should be really easy for a
compute shader to handle- though the more efficiency we can squeeze out of it, the better. The
mirror-style boundary conditions I previously drafted for panda3d-ssbo would probably not
run on many laptops or older/modest home computers.
Force and Velocity Spatial Hashing ODE in RK4 Solver
- FUNCTION ode:
- detect nearby boids using space hash
- calculate force += cohesion - separation
- calculate acc = force / mass (skip this by making boids have constant mass of 1)
- vel += acc
- pos += vel
- return new vel and pos
- detect nearby boids using space hash
START:
- [load prev pos and vel]
- rk4:
- // calculates result of motion across timestep by checking force and recalculating vel across
- // the timestep using the RK4 weighted average method of calculating rates of change.
- k1p, k1v = ode(pos, vel)
- k2p, k2v = ode(pos+k1p/2,vel+k1v/2)
- k3p, k3v = ode(pos+newpos/2,vel+newvel/2)
- k4p, k4v = ode(pos+newpos,vel+newvel)
- vel = (k1v + 2*k2v + 2*k3v + k4v)/6
- pos += vel
- [update new pos and vel]
--- OR ---
Spatial Hashed Force Solver into RK4 Velocity and Position Integrator
- FUNCTION ode(pos, vel, acc):
- vel += acc
- pos += vel
- return new pos and vel
START:
- [load prev pos and acc]
- detect nearby boids using space hash
- calculate force = cohesion - separation
- calculate new acc = force / mass (skip this by making boids have constant mass of 1)
- calculate change in acc = new acc - prev acc
- rk4:
- // now calculates result of motion across timestep using start and end acc values
- k1p, k1v = ode(prev pos, prev vel, prev acc)
- k2p, k2v = ode(prev pos + (k1p - prev pos)/2, prev vel + (k1v - prev vel)/2, prev acc + change in acc / 2)
- k3p, k3v = ode(prev pos + (k2p - prev pos)/2, prev vel + (k2v - prev vel)/2, prev acc + change in acc / 2)
- k4p, k4v = ode(k3p, k3v, new acc)
- new vel = (k1v + 2*k2v + 2*k3v + k4v)/6
- new pos = (k1p + 2*k2p + 2*k3p + k4p)/6
- [update new pos and vel]