A glimpse into a world of 2D particles

This post is inspired by this beautiful youtube video from Tom Mohr (highly recommended watch). It does not really follow the modeling from Tom's video precisely, but it is quite similar. In the process of writing this post I learnt a lot! it felt like it'd been a long time since I did something this fun, therefore I am continuing this project hopefully in a couple of future parts (there is at least going to be a part 2). This series is all about efficiency, modeling, and performance when it comes to web based visualisations. This post is basically my attempt to share what I have learnt along the way and to also document it for myself.
Here's a table of contents for easy scroll to each section What are we looking for? How do we model it? How do we simulate it? (basic implementation) How does it look? (an interactive demo) Tips for better performance Example visualisations What's next? (Part 2)

# What’s the challenge here?

Picture this, you have a bunch of particles (tiny little circles) split into multiple colors. Each color represents a unique type, and the way these particles interact is governed (for the most part) by a fixed matrix, a kind of table where each cell tells you how strongly one type pulls (attraction) or pushes another (repulsion). The rows and columns correspond to the colors, and the number in each cell is how strongly a pair of particles with those colors pull or push one another. don't worry if it all sounds pretty abstract for now, I am going to explain all of this in more details very shortly. Our main goal is to bring this system to life visually and see how these particles interact with different configurations. In a second part to this post coming some time later I will dig deeper into how we can do this in scale with lots of particles and more efficientely in terms of performance and user experience as well.

In a nutshell our objective is to

Model how particles move in every frame of time mathematically
Design an algorithm that implements our modeling frame by frame
Target a frame budget of 16 milliseconds per frame to maintain roughly a stable 60fps performance on Mac Mini M4 Pro hardware when simulating 1000 particles

For the most of this post and this part in particular we will focus on the first two objectives. towards the end of this post we cover the third objective and go even beyond our initial target and aim to make our demo efficient enough for a 60fps performance on a much less powerful hardware such an iPhone 13 Pro, or a Google Pixel 7a.

Modeling the Problem

So far we know we want particles that can move around, they have different colors, each color can attract or repel another color, and that's about it. It is entirely up to us how to mode and then implement these requirements. some modeling choices make it easier for us not only to communicate our final implementation with others (e.g. by use of familiar notations), but to even optimise it further down the road. for instance if we can model our problem in matrix notations, we could potentially use results from Linear Algebra to reduce number of computations in our simulation, or we can utilise GPUs to parallelise some computations. Anyways, in this section we clarify our modeling choices.

Modeling Motion

According to Newton's second law of motion, force is related to acceleration by

Let's assume all particles have unit mass ( ), then for each particle this simplifies to

The acceleration affects the velocity and position of particle over time. a standard and simple way to model these updates is to use the Semi-implicit Euler method at every step, where is the simulation time step. The update equations are:

This models how forces between particles drive their motion by connecting them to acceleration, velocity, and position updates.

In practice, we further apply a small damping factor to the velocity at each time step (pre-computed per particle during initialisation)

where acts as a simple model of energy loss (e.g., from friction or air resistance) and helps maintain numerical stability. In each update step, acceleration is recomputed based on forces, but velocity is accumulated over time. Without damping, even small persistent forces can cause velocity to grow unchecked due to this accumulation. Applying scales down the updated velocity slightly each step, counteracting unrealistic buildup of motion and simulating gradual energy dissipation.

For more on modeling motion I highly recommend this amazing tutorial by Kenichi Yoneda. we can try and build upon this basic motion model later on but for now this works just fine.

Next we look into how we can model pairwise particle forces such that we get attraction and repulsion behaviour based on particle colors.

Modeling Forces

Earlier we briefly hinted at a table to represent pairwise forces. here we make it more explicit by considering it as a matrix called the force matrix. following is an example of this matrix for a system with three particle types

	Red	Green	Blue
Red	0.5	0.3	-0.3
Green	-0.5	0.2	0.5
Blue	0.2	-0.2	0.3

Let's call this matrix :

Each entry represents the force exerted by a particle of color associated with row towards a particle of color associated with column . Diagonal elements (e.g. ) define self-interaction between particles of the same color. For instance in the above example

means red particles attract each other by a force with a magnitude of
means red particles are attracted towards green particles by a force with a magnitude of . (a little weaker than the force from before)
means blue particles repel green ones by a force with a magnitude of

So far for a pair of particles we know how much force is determined by their types, but if this is the only source of force and assuming we have only 2 types that attract each other, after a while all particles collapse into each other. basically we need to model a global repulsive force as well so when particles are too close to one another they repel each other regardless of their type. there are so many possibilities but here we use a very simple one. we model a multiplier which is a function of the distance between particles. let's call it . we use a simple piecewise linear function with a range of where for particles closer than a short distance it gives a negative value regardless of their type (indicating repulsion). In modeling we have 3 key distances

Min Cutoff distance
For particles any closer we clamp the repulsion force to -1
Zero Force distance
A distance at which particles do not apply any force on one another, but any closer they repel and any farther they attract each other
Max Cutoff distance
For particles any farther we assume the force between them is negligible i.e. 0

Here's an illustration of this multiplier function

Force multiplier function

models a source of global repulsion (red part), a limit on how far particles can influence each other, and a simple way to control the distance between neighboring particles

Using , , and for the three distances introduced above, the multiplier in explicit form is:

With and as the row and column of for the colors of particles and , the force strength from particle towards particle is:

Here is the distance between particles and at time , and the sign of determines attraction () or repulsion ().

We can then define force vector , force from particle towards particle as the following

The net force on a particle is the sum of all such vectors i.e. forces coming from all the other particles towards

To understand the vector math above a little better imagine we have only 4 particles , , , and . we want to compute the overal force on particle at a any given time depending on where all the particles are placed and what kinds of attraction and repulsion forces they apply to each other. let's further assume at a given time particle is somewhere in the middle of the others as illustrated in the figure below

2D representation of 4 particles of different colors at a given time

Here

Arrows with solid lines represent (difference) vectors between each pair of particles that include particle . (these are the only pairs that affect )
The arrow with dotted line represents the negated vector between particles and , or more precisely . as you can see in the formula to compute this points towards the direction when or when . since we are are talking about forces applied to particle , this represents a repel force from particle
Dashed arrows represent vectors and (to reduce visual clutter, arrow notation denoting vector quantities has been omitted from the graphics)

Let's further assume green and red both attract orange, but blue repels it (each can have a different constant factor). with these assumptions we can compute , , and and their sum looks like the orange vector below. Plugging in this force into our modeling of motion we can compute in which direction and by how much particle moves and will be positioned at time

Red and green particles both attract orange (green likey has a stronger attraction constant), and blue repels orange. white arrows show all forces applied to particle

centered at

which is its position at time

. orange arrow illustrates the net force.

Boundary Conditions

An important aspect of such a simulation is designing how particles behave at the edges of the 2D world. there are multiple options. Here are three well-known boundary options

Reflecting
Particles bounce off the edges of the boundary
Fixed
Particles are stopped or pinned when they reach an edge of the boundary
Periodic (Wrap-around/toroidal/wrapping)
Particles leaving one side instantly reappear on the opposite side, creating a continuous loop

Each choice has implications for the resulting dynamics and can be selected based on our desired simulation behavior.

Implementation

To better understand how the theoretical model translates into a working simulation, we go one step closer towards the final visualisation by mapping out a flowchart of the simulation logic. this simulation is essentially governed by the modeling we already discussed. so for instance "Compute net force" step is a short way of saying we compute

Interactive Demo

In the following interactive demo you can see how this modeling and simulation logic all comes together towards realising our first two objectives. try playing around with different simulation parameters such as modifying entries of the color matrix, timestep, zero force distance, ... to see in real-time how these particles interact.

Thoughts on Performance

So far we saw the modeling and a demo, but didn't really discuss implementation details. the flowchart from earlier does give us an overal picture however. we have a long running loop and we want each iteration not to take any longer than 16 milliseconds or else it means our framerate drops to below 60fps. part of this has to do with the complexity of the problem itself i.e. the number of particles involved. no matter what we do there is always a large enough number of particles our implementation cannot handle efficiently. so let's focus on our 3rd objective (1k particles at 60fps on a modern Mac Mini) and see what we can do to support simulating that many particles.

Why 60 FPS and scale matter?

Our goal isn’t just to visualize particles, it’s to create a smooth 60fps simulation where you can watch complex behaviors emerge in real time. Think of it like a video game. If the framerate stutters, interactivity almost disappears. Every frame requires calculating forces between every pair of particles, which can quickly tank performance. So basically with each iteration of our implementation we are interested in knowing how many particles can we simulate before framerate drops below 60fps.

Let's briefly discuss the performance of a naive implementation of our modeling. this part requires some familiarity with typescript and p5js. (I am using p5js in instance mode so calls to p5 specific API look like p.createVector instead of createVector

In the spirit of OOP let's start with defining a couple of interfaces that we are going to use in our main code

    
interface Boundary {
    left: number;
    right: number;
    top: number;
    bottom: number;
}

interface PhysicsState {
    position: { 
        x: number,
        y: number
    };
    velocity: {
        x: number,
        y: number
    };
}

interface Force {
    x: number;
    y: number;
}

interface Particle {
    x: number;
    y: number;
    vx: number;
    vy: number;
    type: number;
    velocityDampingMultiplier: number;
}

interface SimulationContext {
    // this many particles are split into three types of colors (almost equally) 
    numParticles: number;
    // radius of each particle
    particleRadius: number;
    // if two particles are closer than this we assume this is their distance
    cutoffMinDistance: number;
    // two particles farther than this do not apply any forces to one another
    cutoffMaxDistance: number; 
    // distance at which particle pairs apply 0 force. any closer they repel each other regardless of their color and any farther they attract/repel each other according to their color.
    zeroForceDistance: number;
    // in seconds, used when computing velocity and position from applied force (F=ma)
    simulationTimestep: number;
    // matrix containing forces between particles based on color
    colorForceMatrix: number[][];
    // simulation boundary (2D box)
    boundary: Boundary;
}

We can then define our update logic like the following

    
function updateParticles(ctxt: SimulationContext) {
    // 1. Computing net force on each particle
    let forcesX = new Array(particles.length).fill(0);
    let forcesY = new Array(particles.length).fill(0);
    for (let i = 0; i < particles.length; i++) {
        let pi = particles[i];
        for (let j = 0; j < particles.length; j++) {
            if (i === j) continue;
            let force = computeForceBetweenParticles(pi, particles[j], ctxt);
            forcesX[i] += force.x;
            forcesY[i] += force.y;
        }
    }
    // 2. Applying motion model based on computed net forces
    for (let i = 0; i < particles.length; i++) {
        let pi = particles[i];
        const currentState: PhysicsState = {
            position: { 
                x: pi.x, 
                y: pi.y 
            },
            velocity: { 
                x: pi.vx,
                y: pi.vy
            }
        };
        const appliedForce: Force = {
            x: forcesX[i],
            y: forcesY[i]
        };
        // Calculate next state (stat at t + 1)
        const nextState = updatePhysics(currentState, appliedForce, pi.velocityDampingMultiplier, ctxt);
        pi.x = nextState.position.x;
        pi.y = nextState.position.y;
        pi.vx = nextState.velocity.x;
        pi.vy = nextState.velocity.y;
        // Applying reflecting boundary conditions (mutates pi position and velocity if needed)
        handleBoundaryConditions(pi, ctxt);
    }
}

The three methods used within updateParticles are as follows

computeForceBetweenParticles
This is where we compute force on particle from particle
updatePhysics
This is where we compute new velocity and position of a particle given a net force vector according to our modeling of motion
handleBoundaryConditions
This is where we constrain particle motion based on reflecting boundary conditions

Here's how they can be implemented

    
function computeForceBetweenParticles(pi: Particle, pj: Particle, ctxt: SimulationContext) {
    let force = p.createVector(pj.x, pj.y).sub(p.createVector(pi.x, pi.y));
    let distance = force.mag();
    // treat any distance less than cutoffMinDistance as cutoffMinDistance
    if (distance < ctxt.cutoffMinDistance) distance = ctxt.cutoffMinDistance;
    // applying our modeling for attraction/repulsion
    let strength = 1;
    let multiplier = 1;
    if (distance <= ctxt.zeroForceDistance) {
        // cutoffMinDistance <= distance <= zeroForceDistance 
        // distance too close -> global repulsion to prevent collision
        multiplier = 1 - p.map(distance, ctxt.cutoffMinDistance, ctxt.zeroForceDistance, 0, 1);
        // fixed repulsion force with magnitude of 1 (regardless of particle colors)
        strength = multiplier * -1;
    } else if (distance <= ctxt.cutoffMaxDistance) {
        // zeroForceDistance < distance <= cutoffMaxDistance 
        multiplier = p.map(distance, ctxt.zeroForceDistance, ctxt.cutoffMaxDistance, 0, 1);
        strength = multiplier * ctxt.colorForceMatrix[pi.type][pj.type];
    } else {
        // distance > cutoffMaxDistance so force between the pair is considered negligible
        strength = 0;
    }
    return force.setMag(strength);
}

Then we have

    
const updatePhysics = (
    current: PhysicsState,
    force: Force,
    particleVelocityDampingMultiplier: number,
    ctxt: SimulationContext
): PhysicsState => {
    const dt = ctxt.simulationTimestep;
    // 1. Calculate acceleration components (assuming unit mass i.e. m=1 so a=F/m=F)
    const ax = force.x;
    const ay = force.y;
    // 2. Update velocity components (v = v0 + a*dt)
    let vxNew = current.velocity.x + ax * dt;
    let vyNew = current.velocity.y + ay * dt;
    // 3. Apply damping to updated velocities to avoid runaway speeds
    vxNew *= particleVelocityDampingMultiplier;
    vyNew *= particleVelocityDampingMultiplier;
    // 4. Update position components (x = x0 + v*dt)
    const xNew = current.position.x
        + vxNew * dt;
    const yNew = current.position.y
        + vyNew * dt;
    return {
        position: { 
            x: xNew,
            y: yNew
        },
        velocity: {
            x: vxNew,
            y: vyNew,
        }
    };
}

There are several options how to define boundary conditions and each has its own consequences and meanings. In this post for simplicity we only consider reflecting boundary conditions (without any energy loss) which can be implemented simply as below

    
function handleBoundaryConditions(particle: Particle, ctxt: SimulationContext) {
    // bouncing off the walls
    if (particle.x > ctxt.boundary.right - ctxt.particleRadius) particle.vx *= -1;
    if (particle.x < ctxt.boundary.left + ctxt.particleRadius) particle.vx *= -1;
    if (particle.y > ctxt.boundary.bottom - ctxt.particleRadius) particle.vy *= -1;
    if (particle.y < ctxt.boundary.top + ctxt.particleRadius) particle.vy *= -1;
    // making sure particles stay within the boundary
    particle.x = p.constrain(particle.x, ctxt.boundary.left + ctxt.particleRadius, ctxt.boundary.right - ctxt.particleRadius);
    particle.y = p.constrain(particle.y, ctxt.boundary.top + ctxt.particleRadius, ctxt.boundary.bottom - ctxt.particleRadius);
}

In part 2 we will also explore implementation details of periodic boundary conditions.

Now that we have familiarised ourselves with the code, first point that comes to mind is that within updateParticles we are looking at every pair twice. once computing force from particle towards and vice versa. these forces are not necessarily equal as you'd expect because the colorForceMatrix is not necessarily a symmetric one. However, one major (heavy) computation here is computing the distance between the two particles and in computing both and this distance is the same. So we better only compute it once per pair. Let's refactor updateParticles based on this observation (notice the highlighted lines)

    
    // ...
    for (let i = 0; i < particles.length - 1; i++) {
        for (let j = i + 1; j < particles.length; j++) {
            // we compute not just the force from j to i but also the one from i to j 
            // (this way we compute distance between i and j only once)
            let force_ij = p.createVector(particles[j].x, particles[j].y).sub(p.createVector(particles[i].x, particles[i].y));
            let force_ji = force_ij.copy();
            let distance = force_ij.mag();
            if (distance < ctxt.cutoffMinDistance) distance = ctxt.cutoffMinDistance;
            let strength_ij = 1;
            let strength_ji = 1;
            let multiplier = 1;
            if (distance <= ctxt.zeroForceDistance) {
                multiplier = 1 - p.map(distance, ctxt.cutoffMinDistance, ctxt.zeroForceDistance, 0, 1);
                strength_ij = -1 * multiplier; // we consider max repulsion to be -1
                strength_ji = strength_ij;
            } else if (distance <= ctxt.cutoffMaxDistance) {
                multiplier = p.map(distance, ctxt.zeroForceDistance, ctxt.cutoffMaxDistance, 0, 1);
                strength_ij = multiplier * ctxt.colorForceMatrix[particles[i].type][particles[j].type];
                strength_ji = multiplier * ctxt.colorForceMatrix[particles[j].type][particles[i].type];
            } else {
                strength_ij = 0;
                strength_ji = 0;
            }
            force_ij.setMag(strength_ij);
            force_ji.setMag(strength_ji);
            forcesX[i] += force_ij.x;
            forcesY[i] += force_ij.y;
            forcesX[j] -= force_ji.x; // opposite force direction
            forcesY[j] -= force_ji.y;
        }
    }
    // continue with applying motion model using forcesX[i] and forcesY[i] for every particle
    // ...

With this simple change we reduce redundant force calculations by roughly half. This is just the start.

Now let's look at something a little more subtle. within updateParticles we use p.createVector and a couple of related methods such as .sub(), .mag(), and .setMag(). This essentialy incurs some object creation overhead and adds pressure to javascript engine's garbage collection. We could simply get rid of these and use direct numerical operations and primitive types. let's see how that looks

    
    // ...
    for (let i = 0; i < particles.length - 1; i++) {
        for (let j = i + 1; j < particles.length; j++) {
            const dx = particles[j].x - particles[i].x;
            const dy = particles[j].y - particles[i].y;
            let distance = Math.sqrt(dx * dx + dy * dy);
            if (distance < ctxt.cutoffMinDistance) distance = ctxt.cutoffMinDistance;
            let strength_ij = 1;
            let strength_ji = 1;
            let multiplier = 1;
            if (distance <= ctxt.zeroForceDistance) {
                multiplier = 1 - p.map(distance, ctxt.cutoffMinDistance, ctxt.zeroForceDistance, 0, 1);
                strength_ij = -1 * multiplier; // we consider max repulsion to be -1
                strength_ji = strength_ij;
            } else if (distance <= ctxt.cutoffMaxDistance) {
                multiplier = p.map(distance, ctxt.zeroForceDistance, ctxt.cutoffMaxDistance, 0, 1);
                strength_ij = multiplier * ctxt.colorForceMatrix[particles[i].type][particles[j].type];
                strength_ji = multiplier * ctxt.colorForceMatrix[particles[j].type][particles[i].type];
            } else {
                strength_ij = 0;
                strength_ji = 0;
            }
            forcesX[i] += (dx / distance) * strength_ij;
            forcesY[i] += (dy / distance) * strength_ij;
            forcesX[j] -= (dx / distance) * strength_ji;
            forcesY[j] -= (dy / distance) * strength_ji;
        }
    }
    // continue with applying motion model using forcesX[i] and forcesY[i] for every particle
    // ...

Even though in the first step we directly targetted redundant computations, performance effects from this new update were actually more noticable. This shows how much time savings within nested loops can impact performance.

The next small but useful step (V3) is to hoist simulation constants (cutoff distances, the force matrix) out of the nested loop into local const variables before the loop starts. Accessing a local variable is cheaper than a property lookup inside a hot loop that runs millions of times per frame.

V4 then goes one step further and looks at how we store particle data itself. OOP is nice and all but to speed up further we can look for more memory lookup efficiency. This can be achieved by storing particle data in a flat array instead of an array of objects. Here's how it looks

    
    const PARTICLE_INFO_SIZE = 6; // x, y, vx, vy, type, velocityDampingMultiplier
    const X = 0, Y = 1, VX = 2, VY = 3, TYPE = 4, DAMPING_FACTOR = 5;
    let particles = new Float32Array(ctxt.numParticles * PARTICLE_INFO_SIZE);

    function updateParticles(ctxt: SimulationContext) {
        // 0. we precompute frequently used constants within the loops
        const minDist = ctxt.cutoffMinDistance;
        const maxDist = ctxt.cutoffMaxDistance;
        const zeroForceDist = ctxt.zeroForceDistance;
        const colorForceMatrix = ctxt.colorForceMatrix;
        const n = ctxt.numParticles;
        // ...
        for (let i = 0; i < n - 1; i++) {
            let baseI = i * PARTICLE_INFO_SIZE;
            let xi = particles[baseI + X];
            let yi = particles[baseI + Y];
            let typeI = particles[baseI + TYPE];
            for (let j = i + 1; j < n; j++) {
                let baseJ = j * PARTICLE_INFO_SIZE;
                let xj = particles[baseJ + X];
                let yj = particles[baseJ + Y];
                let typeJ = particles[baseJ + TYPE];
                const dx = xj - xi;
                const dy = yj - yi;
                let distance = Math.sqrt(dx * dx + dy * dy);
                if (distance < minDist) distance = minDist;
                let strength_ij = 1;
                let strength_ji = 1;
                let multiplier = 1;
                if (distance <= zeroForceDist) {
                    multiplier = 1 - p.map(distance, minDist, zeroForceDist, 0, 1);
                    strength_ij = -1 * multiplier; // we consider max repulsion to be -1
                    strength_ji = strength_ij;
                } else if (distance <= maxDist) {
                    multiplier = p.map(distance, zeroForceDist, maxDist, 0, 1);
                    strength_ij = multiplier * colorForceMatrix[typeI][typeJ];
                    strength_ji = multiplier * colorForceMatrix[typeJ][typeI];
                } else {
                    strength_ij = 0;
                    strength_ji = 0;
                }
                forcesX[i] += (dx / distance) * strength_ij;
                forcesY[i] += (dy / distance) * strength_ij;
                forcesX[j] -= (dx / distance) * strength_ji;
                forcesY[j] -= (dy / distance) * strength_ji;
            }
        }
        // continue with applying motion model using forcesX[i] and forcesY[i] for every particle
        // ...
    }

Other methods of course need to be adapted as we no longer have objects of type Particle.

With these optimisations we can easily render slightly over 2000 particles at 60fps on a Mac Mini M4 Pro which was our initial objective.

Examples

In this section I want to share a couple of examples of how this simulation can be used to create interesting visual effects. in all of these examples the simulation parameters are kept the same except for the color force matrix. Initially particles are distrubuted around a circle and the boundary conditions make it so they are constrained within a white rectangular box.

A particle simulation where the color force matrix (

) is set to zero. therefore no pulling or pushing forces exist between particles due to their colors.

A simulation where

is set an identity matrix. therefore particles of the same color attract each other while particles of different colors do not interact.

A simulation where all entries of

are set to 1. therefore particles of all colors attract each other equally.

A simulation where all entries of

are initially set to 1 and then diagonal entries are inverted. therefore particles of the same color repel each other while particles of different colors attract each other.

What's coming in part 2?

In part 2 we break the barrier with spatial partitioning and move physics off the main thread using Web Workers. A WebGPU compute path pushes the simulation further, reaching 10,000 particles at 60fps on mobile. Future parts may also explore other directions such as WebAssembly or matrix-based formulations, so keep an eye out for more updates!