Chaotic Systems, Realtime Particle Simulations and 3D Rendering

Simulating multiple independent particles following a chaotic ODE, and rendering particle trails with Three.js.

Contents

Part 1: Realtime Simulator
Part 2: Trail Rendering
Final Result

Part 1: Realtime Simulator

Dynamical Systems

Imagine a set of particles moving around in 3-dimensional space, each following the dynamics of an ODE:

\dot{\bm{x}} = f(t, \bm{x})

The behavior of each particle is then determined by the form of the dynamics. And there are of course a lot of different possibilities. A common way to look at the behavior is to look at the various equilibrium points and how they make the system behave.

A linear autonomous system given as:

\dot{\bm{x}} = \bm{A} \bm{x}

for some square matrix $\bm{A}$ , will have its global behavior determined by the type of equilibrium point at the origin. But basically, most trajectories will either converge to the origin if the linear system is stable, or diverge to infinity for an unstable system (marginal stability being the third case).

For visualization purposes, neither of the above are particularly interesting to look at. Which is where nonlinear dynamical systems come in, giving more interesting visuals. Global convergence to a single point, a periodic orbit, or divergence is still possible, which again depends on the dynamics.

Chaotic Systems

For a dynamical system to be chaotic in some subspace it's enough to satisfy these requirements:

Topological transitivity
Dense periodic orbits

which is not necessarily trivial to show for a general system. For a more informal explanation, the first requirement kind of means that trajectories will come arbitrarily close to visiting every part of the subspace. And the second requirement means that most trajectories are always very close to an unstable periodic orbit, which provides some structure to the way they move around.

Robert Devaney's original formulation included one more requirement, which is what many people often associate with chaotic systems: sensitive dependence on initial conditions. This approximately means that two trajectories that start very close together will eventually start diverging from each other at an exponential rate (measured by what is called a Lyapunov exponent). It was later shown by Banks et al. that this requirement can be proved to be satisfied if the first two requirements are fulfilled, and is therefore unnecessary to include in the definition.

And in even simpler words, chaotic systems are cool to visualize because all trajectories just move around endlessly, but still contain some structure. This subspace is called an attractor if trajectories outside converge into it. And many chaotic systems result in attractors containing a fractal dimension, a little because the trajectories almost work like space-filling curves, which have the interesting name of a strange attractor.

The Lorenz System

Chaotic systems can be both discrete time, and continuous time. A lot of textbooks like to start with talking about discrete-time systems to introduce chaos, but I am personally more fascinated by continuous time systems. The Lorenz system was historically one of the first examples that explicitly demonstrated chaotic behavior, and the concept of sensitive dependence on initial conditions.

The system is defined as the following 3-dimensional nonlinear ODE:

\begin{align*} \dot{x} &= \sigma (y - x) \\ \dot{y} &= x (\rho - z) - y \\ \dot{z} &= x y - \beta z \end{align*}

with three parameters: $\sigma, \rho, \beta$ . The specific values of these parameters will greatly influence the behavior of the overall system, containing multiple bifurcation points. The most common parameters to use for chaotic behavior are: $\sigma = 10, \rho = 28, \beta = 8 / 3$ , which gives three unstable equilibrium points: the origin is always an equilibrium point, and two others in the centers of the two discs of the attractor.

Lorenz attractor — Lorenz butterfly attractor. Name is probably related partially to the looks, partially for the butterfly effect.

Numerically Integrating Chaotic Systems

One thing that makes chaotic systems interesting from a more applied math perspective is that they are very hard to solve numerically. Computers can't represent real numbers, and rely on approximating them, typically with floating point numbers. But this leads to some inaccuracies from the quantization, and one of the key properties of chaotic systems is that a small error compounds to an exponentially increasing difference. This means that any numerical solution will eventually be too inaccurate to be of use after some integration time.

So this could make it seem like fancy integration methods don't work, but they can still give a better prediction horizon than a more naive integrator. To see this, look at how error accumulates in chaotic systems. Start with trajectories $\bm{x}_1(t)$ and $\bm{x}_2(t)$ such that:

\bm{\epsilon}(t) = \bm{x}_1(t) - \bm{x}_2(t)

is the difference between them. If we assume that the initial error $\bm{\epsilon}(0) = \bm{x}_1(0) - \bm{x}_2(0)$ is small, then the Lyapunov exponent of the system $\lambda$ says that:

|| \bm{\epsilon}(t) || \approx || \bm{\epsilon}_0 || e^{\lambda t}

We can then define a maximum acceptable error: $\epsilon_{\textrm{max}}$ , and manipulate the equation above into this:

t_{\textrm{max}} = \frac{1}{\lambda} \ln{\left( \frac{\epsilon_{\textrm{max}}}{|| \bm{\epsilon}_0 ||} \right)}

so if we can reduce $\epsilon_{\textrm{max}}$ from fancy integration methods, it will contribute to extending the prediction horizon $t_{\textrm{max}}$ . Although it is logarithmic in the error, so increasing the horizon requires exponentially more effort from the integrator.

Also note that integration accuracy doesn't actually matter for visualizing it here, as all trajectories stay on the attractor anyways. Or put another way, it does matter for an individual trajectory, but not really for visualizing the attractor itself with a set of trajectories, which is kind of what this is building up to.

Diverging trajectories — Top plot shows how two trajectories with almost the same initial conditions start diverging after some time. Bottom plot shows the norm of the difference between them. Notice how the difference has a section where it increases linearly on a log scale.

Part 2: Trail Rendering

Timestep Interpolation

To make things look more interesting visually, each particle is rendered with a trail. To do this, save not just the particle position, but also the history. This is done in an efficient way with ring buffers and lots of index manipulations to avoid having to allocate memory dynamically at runtime.

So every iteration of the particle simulation, update the current particle index and insert the position from the simulator. It currently uses the Dormand-Prince 5 method as an adaptive timestep integration algorithm, but the simulator is still running on a fixed timestep in between particle position saves.

To draw in the smoothest way possible, the rendering should use the same update frequency as the refresh rate of the current display. Because of the fixed physics frequency, these will rarely align. So what the overall system does is keep track of a timestep accumulator that increments with the deltatime of the rendering loop. Whenever the accumulator is greater than a physics timestep, then do an iteration of the simulation state. Depending on the display refresh rate and the chosen fixed physics frequency, it can either do multiple physics updates per rendered frame, or multiple renders per each physics update. On average.

Because the timings don't work out exactly, the previous simulation state is also saved, and the exact state to be rendered is found by linearly interpolating between those two states.

Different frequencies — Blue dots represent the physics updates. Orange dots represent the rendering updates. So for a given orange dot, interpolate the physics state between the two physics states around it.

Drawing Lines

Turns out to be non-trivial.

Using simple line drawing algorithms like Bresenham kind of gives lines with a width of 1. Which is not really what this rendering is trying to do. There are different ways to add some width, but they all involve constructing a mesh. For example a series of rectangles like what is done here. This is also called a ribbon.

For each two points in a particle trail, the goal is then to draw a rectangle. This is done by finding the delta position vector, finding the camera position vector and taking the cross product of these two. This new vector points in a direction which is perpendicular to both the camera and to the direction of the particle.

Given the current and the previous particle positions, adding and subtracting this new perpendicular vector then gives four points each making up the corner of the rectangle. It is also easy to scale the width like this.

Shaders

Computing the ribbon mesh for every particle, at every iteration is very slow. It is possible to construct the ribbon mesh with a ring buffer as well, and only update the latest. But we are drawing on the GPU anyways, and GPUs work best when the amount of data transfers are minimized. So it is also possible to construct the ribbon mesh directly on the GPU using a vertex shader.

A vertex shader is kind of meant to be run for every point in a mesh, and transform the position so it makes sense with respect to the camera. But they can be used in more creative ways. Here the particle positions are sent to VRAM by embedding them in a texture. And the points in the mesh are instead based on particle and trail indexes, along with which side of the rectangle they are located. This means that the vertex shader runs once for every combination of particle index, trail index and side of the rectangle.

Colors

As the colors of each particle trail is based on the positions and velocities, they are computed in the vertex shader and passed along to the fragment shader. There's also a lot of additional color state variables just meant to make it look more interesting that are passed in as uniforms.

The two main ways that colors are set is that the angle between the x and y components of the particle velocity is mapped to a hue value, and that the alpha value is based on the age of the position so each trail fade out. The hue value is then stretched to cover some range of 0 to 360 as with HSL colors, while the saturation and light are based on a global state that changes randomly. These are then mapped to RGB in the shader.

Final Result

Orbit controls are also implemented, so try dragging it around to see from other angles, or zoom in and out.

The full source code can be found here implemented with TypeScript and GLSL.