Intro
These aren’t comprehensive lecture notes or a textbook. They’re simply my attempt to organize everything I have to learn to understand fluid simulation beyond just copying code.
What started as “I’ll build a simple fluid simulator this weekend” quickly turned into weeks of relearning calculus, numerical methods, and linear algebra before I could even understand the papers I was trying to read.
I’m still very much learning. As I continue exploring the topic, I’ll keep adding notes on the mathematics, numerical methods, and simulation techniques I find useful. If you’re trying to learn this stuff too, maybe they’ll save you a few of the detours I had to take.
Real-Time Fluid Dynamics for Games
Before writing Real-Time Fluid Dynamics for Games, Jos Stam published Stable Fluids in 1999, introducing the ideas that made his approach to fluid simulation so influential. From what I understand so far, this later paper is a much smaller and more approachable introduction to those same ideas, aimed at getting something running quickly rather than diving into all the theory.
The goal isn’t scientific accuracy. Instead, it’s to produce visually convincing smoke and fluid motion while remaining stable and fast enough for real-time applications. Better yet, the paper includes the complete C implementation, making it possible to follow both the mathematics and the code side by side.
That’s why I’m starting here. Rather than treating the paper as something to blindly implement, I’m using it as an excuse to finally understand the mathematics behind fluid simulation. I’ll almost certainly branch into other papers, books, and techniques as I learn more, but this feels like a good place to start.
Navier-Stokes equations
The big bad wolf. The million-dollar equation (literally — there’s a Millennium Prize for anyone who can solve the Navier–Stokes existence and smoothness problem).
Stam writes, in the most compact vector form, the equation for the velocity :
$$ \frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} $$
Here $\mathbf{u}$ is the velocity field — a vector field telling you the velocity of the fluid at each point in space. The way I see it, this equation is kind of saying: the change in velocity over time (the left side) comes from three main things: how the fluid moves itself around, how it diffuses or spreads out due to viscosity ($\nu \nabla^2 \mathbf{u}$ term), and whatever external forces are acting on it ($\mathbf{f}$).
The second equation describes how a density field (like a continuous smoke density cause we are not modeling every particle) moves through the velocity field.
$$ \frac{\partial \rho}{\partial t} = -(\mathbf{u} \cdot \nabla)\rho + \kappa \nabla^2 \rho + S $$
The main idea is that the density $\rho$ gets carried around by the velocity field (first term), it diffuses with some diffusion constant $\kappa$, and there might be sources $S$ that add more density (like smoke being emitted).
In Fluid Simulation for Computer Graphics, Bridson phrases the Navier-Stokes equations like so:
“Most fluid flow of interest in animation is governed by the famous in-compressible Navier-Stokes equations, a set of partial differential equations that are supposed to hold throughout the fluid”
Momentum Equation:
$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho} \nabla p + \mathbf{g} + \nu \nabla^2 \mathbf{u} $$
- $\mathbf{u}$ — velocity of the fluid
- $p$ — pressure (force per unit area that the fluid exerts on anything)
- $\rho$ — density of the fluid (mass per unit volume $kg/m^3$)
- $\mathbf{g}$ — gravity and any additional external forces
- $\nu$ — kinematic viscosity (how “thick” or resistant to stirring the fluid is)
Incompressibility Condition:
$$ \nabla \cdot \mathbf{u} = 0 $$
This expresses that the fluid is incompressible — its volume doesn’t change, and there’s no net gain or loss of fluid at any point.
So, depending on what assumptions you make — incompressible vs. compressible flow, constant viscosity, ignoring temperature, etc. — you can write these equations in a bunch of different ways. The ones Stam uses are simplified for real-time visual simulation: he’s not trying to simulate airplane aerodynamics, just something that looks like smoke or water moving naturally.
For now, I’m not worrying about all the physics behind them — like where exactly they come from or what all the simplifications mean. I just want to understand how the equations turns into the little C code blocks in Stam’s paper.
Solver Basic Structure
The discrete representation is a 2D grid where each cell stores quantities like velocity and density. Values are assumed to be constant within each cell and represented at the center (later Bridson explains methods where not storing at the center is a much better idea).
Bridson calls this Eulerian Viewpoint TODO continue
We surround our grid simulation region with an aditional layer of cells that act as a buffer region to handle boundary conditions. In practice we can store everything in a one dimensional grid $size = (N+2)*(N+2)$, the paper describes how.
Our solver begins with an initial state for both velocity and density fields, then updates them based on interactions over time. In practice, this means we repeatedly take the current state of the fluid, apply forces, inject new densities, and step the simulation forward by a small time interval $dt$ which in our case will be our frame rate.
Moving Densities
Stam starts by focusing on the simpler of the two equations — the one describing how a density field moves through a fixed velocity field (i.e., assuming the velocity doesn’t change with time).
The density equation:
$$ \frac{\partial \rho}{\partial t} = -(\mathbf{u} \cdot \nabla)\rho + k \nabla^2 \rho + S $$
Each term on the right-hand side contributes differently to how the density evolves:
- The first term means the density follows the velocity field — it gets advected.
- The second term represents diffusion, how density spreads out over time.
- The third term adds sources, meaning new density gets introduced into the system.
Stam’s solver tackles these three effects every time step, but in reverse order:
first adding sources, then diffusing, and finally advecting the density through the velocity field.
The source term is the simplest to handle — for each cell, the new density is increased by the amount added by sources during the time step $dt$; In my implementation I simply set the density at the clicked grid cell to 1.
That’s the first building block — adding density into the world. The interesting behavior comes next, when the density starts to move and spread.
Diffuse Bad
Summary of Stam’s Explanation
In Stam’s paper, the diffusion step accounts for the way density spreads across the grid at a certain rate, diff. When diff > 0, density will naturally distribute among neighboring cells. Each cell exchanges density with its four direct neighbors — left, right, top, and bottom.
For each cell, the net density change can be approximated by:
$$ x_0[IX(i-1,j)] + x_0[IX(i+1,j)] + x_0[IX(i,j-1)] + x_0[IX(i,j+1)] - 4x_0[IX(i,j)] $$
A basic diffusion solver computes this exchange for every cell and adds it to the existing value, resulting in an implementation like this:
void diffuse_bad ( int N, int b, float *x, float *x0, float diff, float dt ) {
int i, j;
float a = dt * diff * N * N;
for ( i = 1 ; i <= N ; i++ ) {
for ( j = 1 ; j <= N ; j++ ) {
x[IX(i,j)] = x0[IX(i,j)] + a * (
x0[IX(i-1,j)] + x0[IX(i+1,j)] +
x0[IX(i,j-1)] + x0[IX(i,j+1)] -
4 * x0[IX(i,j)]
);
}
}
set_bnd(N, b, x);
}
The routine set_bnd() handles the boundary cells, which will be discussed later.
Understanding the Diffusion Equation
If we only account for density diffusion, we obtain the following partial differential equation:
$$ \frac{\partial p}{\partial t} = k \nabla^2 p $$
This is a famous equation — it’s equivalent to the heat equation for a constant diffusion rate $k$.
From Wikipedia: “It describes the macroscopic behavior of many micro-particles in Brownian motion, resulting from the random movements and collisions of the particles.”
It sounds cool, whatever that means.
The Laplacian operator $\nabla^2$ essentially measures how much a point differs from its surrounding neighbors. Personally, I like to think of it as the divergence of the gradient (it can be written $\nabla \cdot \nabla f$).
Given a scalar field of density, the gradient gives us a vector field where the vectors point in the direction of greatest increase — away from low-density areas and toward higher-density ones. The divergence measures how much “outflow” there is at a point:
- Negative divergence → more flow entering than leaving (the region becomes denser).
- Positive divergence → more flow leaving than entering (the region loses density).
- Zero divergence → flow is balanced.
So, the Laplacian, analogous to a second derivative, will be positive where flow is diverging (losing density) and negative where it’s converging (gaining density). It tells us how “peaky” or “valley-like” a point is in the field.
Looking again at our PDE:
$$ \frac{\partial p}{\partial t} = k \nabla^2 p $$
I don’t know how it was originally derived, but intuitively, it says: the rate of change of density over time depends on how density spreads across space. Points with higher density will naturally flow toward areas with lower density.
(There’s a great 3Blue1Brown video that visually explains this beautifully: Diffusion equation)
Solving the Diffusion Equation
What we care about is: given some density $p$ at time step $t$, how do we find the density at $t + \Delta t$?
For this, we use the simplest numerical integrator — the explicit Forward Euler method.
From the definition of the derivative:
$$ \frac{dx}{dt} = \lim_{\Delta t \to 0} \frac{x(t+\Delta t) - x(t)}{\Delta t} $$
We approximate it using a small time step $\Delta t$:
$$ \frac{dx}{dt} \approx \frac{x(t+\Delta t) - x(t)}{\Delta t} $$
Rearranging, we get our update rule:
$$ x(t+\Delta t) = x(t) + \Delta t \frac{dx}{dt} $$
In our case, we don’t directly know $\frac{dp}{dt}$, but our diffusion equation tells us:
$$ \frac{dp}{dt} = k \nabla^2 p $$
Since the Laplacian of the density field, $\nabla^2 p$, is only dependent on the density itself, we just need a way to compute it.
Discretizing the Laplacian
The Laplacian operator can be expanded as:
$$ \nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} $$
We can discretize the second derivative using the central difference formula (see background section on numerical differentiation):
$$ \frac{d^2 f}{dx^2} \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} $$
Substituting into the Laplacian, we get:
$$ \nabla^2 f \approx \frac{f(x+h, y) - 2f(x, y) + f(x-h, y)}{h_x^2} + \frac{f(x, y+h) - 2f(x, y) + f(x, y-h)}{h_y^2} $$
Since we’re using a square grid ($h_x = h_y = h$), this simplifies to:
$$ \nabla^2 f \approx \frac{f(x+h, y) + f(x-h, y) + f(x, y+h) + f(x, y-h) - 4f(x, y)}{h^2} $$
Notice how this matches exactly the expression Stam uses for the diffusion process.
It doesn’t come out of thin air — it’s the result of combining several deep mathematical concepts into something surprisingly simple and elegant.
Putting It All Together
Our update rule for a single cell becomes:
$$ p(t + \Delta t) = p(t) + \Delta t k \frac{p(x+h, y) + p(x-h, y) + p(x, y+h) + p(x, y-h) - 4p(x, y)}{h^2} $$
This looks almost identical to Stam’s pseudo-code.
The only detail to note is that $h$ is the cell size — since each cell has constant density, we must look at neighboring cells to estimate the spatial derivative.
Stam combines $\Delta t$, $k$, and $1/h^2$ into a single constant (the variable a in his code). Because this value doesn’t change during each iteration (or frame), it can be computed once and reused:
$$ a = \Delta t * k * N^2 $$
Here, we assume the grid size is 1, so $h = 1/N$ and $1/h^2 = N^2$.
That’s it — that’s all you need to understand the diffusion step of the simulation. Rejoice!
To experiment with the simulation below, click and drag on the grid to set density values. You’ll find:
- a Play/Pause button,
- a Step button (advances the simulation one step),
- a Refresh button (resets all densities to zero), and
- a slider to control the diffusion coefficient $k$.
Diffuse Good
In the last section we got a minimal diffusion simulation up and running. Pretty cool right?!
But if you tried bumping up the diffusion coefficient, you probably noticed the simulation freaks out — densities start to oscillate, blow up and form weird flickering patterns.
As Stam says in his paper:
“For large diffusion rates the density values start to oscillate, become negative and finally diverge, making the simulation useless. This behavior is a general problem that plagues unstable methods.”
Stam is specifically referring to the fact that the Forward Euler method for numerically integrating the diffusion equation is unstable. What that means is, in the scheme:
$$ x(t + \Delta t) = x(t) + \Delta t \frac{dx(t)}{dt} $$
for a large enough time step $\Delta t$, the values of $x$ will eventually blow up — growing to infinity.
For the diffusion equation $f(x) = k \nabla^2 x$, this instability happens because each update depends entirely on the previous state. Small numerical errors get amplified over time, especially when $\Delta t$ or $k$ (the diffusion rate) is large.
In theory, we could just reduce $\Delta t$ until it stabilizes, but that would mean doing many more iterations per frame — not very practical if we want smooth, real-time behavior.
Actually explaining why the method is unstable is a bit technical, and you can read more about it in the Numerical Integration of ODEs section (TODO write this).
Backward Euler to the rescue
To solve this problem, Stam switches from Forward Euler to Backward Euler, which is an implicit method.
Instead of estimating the future state using the derivative of the current state, Backward Euler uses the derivative of the future state — the one we’re trying to compute.
$$ \frac{dx(t+\Delta t)}{dt} \approx \frac{x(t+\Delta t) - x(t)}{\Delta t} $$
Rearranging:
$$ x(t + \Delta t) = x(t) + \Delta t \frac{dx(t+\Delta t)}{dt} $$
This is what makes the method implicit — the future value $x(t + \Delta t)$ appears on both sides of the equation. We can’t just compute it directly; we have to solve for it.
That might sound like extra work (and it is), but the big advantage is stability: implicit methods like Backward Euler remain stable no matter how large $\Delta t$ is. You can take huge time steps, crank up the diffusion coefficient, and the simulation will still behave.
Let’s see what Backward Euler looks like when applied to the diffusion equation.
The continuous diffusion equation can be written as:
$$ \frac{dx}{dt} = k \nabla^2 x $$
If we apply Backward Euler, we replace the derivative term using the future value $x(t + \Delta t)$:
$$ x(t + \Delta t) = x(t) + \Delta t k \nabla^2 x(t + \Delta t) $$
Rearranging:
$$ x(t) = x(t + \Delta t) - \Delta t k \nabla^2 x(t + \Delta t) $$
And this is exactly what Stam writes in his code (just in discrete form):
x0[IX(i,j)] = x[IX(i,j)] - a * (x[IX(i-1,j)] + x[IX(i+1,j)] + x[IX(i,j-1)] + x[IX(i,j+1)]
- 4 * x[IX(i,j)]
);
Remember $a = \Delta t * k * N^2$ is the scaled diffusion strength.
Because $x(t + \Delta t)$ appears on both sides of the equation — directly and inside the Laplacian — this is an implicit equation. In other words, we can’t just “plug and chug” like before; we have to rearrange and solve for $x(t + \Delta t)$ itself.
Solving the system of equations
What we’re dealing with here is a linear system of equations — we have one equation per cell of our grid.
We can express the problem in matrix form:
$$ A x = b $$
or, more specifically for our diffusion step,
$$ (I - \Delta t k L) x^{n+1} = x^{n} $$
where:
- $b = x^{n}$ is a vector containing our previous density grid values,
- $x = x^{n+1}$ is the vector of unknowns (the new densities we want to compute), and
- $A = (I - \Delta t k L)$ is our matrix of coefficients, with $L$ being the discrete Laplacian matrix.
If we were to solve this explicitly, we’d write:
$$ x = A^{-1} b $$
This is what Stam refers to in his paper:
“This is a linear system for the unknowns $x[IX(i,j)]$. We could build the matrix for this linear system and then call a standard matrix inversion routine. However, this is overkill for this problem because the matrix is very sparse: only very few of its elements are non-zero. Consequently, we can use a simpler iterative technique to invert the matrix.”
Indeed, a quick search shows that most matrix inversion algorithms have a complexity of $O(n^3)$ for an $n \times n$ matrix — which is way too costly for large grids.
Stam instead takes advantage of the structure and sparsity of our coefficient matrix to solve the equations in a much more straightforward and efficient manner.
1D diffusion
For simplicity’s sake let’s suppose instead of a 2D grid we are working with a 1D rod with only 4 cells.
For 1D, the discrete Laplacian is:
$$ \nabla^2 x_i = x_{i-1} + x_{i+1} - 2x_i $$
Substituting this into the Backward Euler equation (to avoid filling the screen with math notation I’ll say $x^{n+1} = x$):
$$ x^n_i = x_i - a(x_{i-1} + x_{i+1} - 2x_i) $$
Expanding terms:
$$ x^n_i = (1 + 2a)x_i - a x_{i-1} - a x_{i+1} $$
In matrix form:
$$ \begin{bmatrix} 1+2a & -a & 0 & 0 \\ -a & 1+2a & -a & 0 \\ 0 & -a & 1+2a & -a \\ 0 & 0 & -a & 1+2a \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} x^n_0 \\ x^n_1 \\ x^n_2 \\ x^n_3 \end{bmatrix} $$
which matches the general form:
$$ (I - \Delta t k L) x^{n+1} = x^{n} $$
At the boundaries, I’m assuming cells outside our domain (like $x_{-1}$ or $x_4$) have value 0. Boundaries are at the edge of my procupations, I’ll leave it for future me to go more in depth on how different boundary conditions affect the problem.
Observing the Matrix
Looking at matrix $A$, we can note the following:
- As Stam said, it’s very sparse — most coefficients are zero (TODO check formally, a sparse matrix is one where the majority of elements are zero).
- It’s diagonally dominant — since $a > 0$, the diagonal entries $(1 + 2a)$ are larger than the sum of the magnitudes of the other entries in each row.
This diagonal dominance is why the Gauss-Seidel relaxation method works here — it’s an iterative technique that approximates the solutions for diagonally dominant systems.
(See the Background: Numerical Solvers for Linear Systems section for more detail.)
Solving with Gauss-Seidel
Let’s rewrite each equation to isolate $x_i$:
$$ \begin{aligned} x_0 &= \frac{x^n_0 + ax_1 + 0x_2 + 0x_3}{1 + 2a} \\ x_1 &= \frac{x^n_1 + ax_0 + ax_2 + 0x_3}{1 + 2a} \\ x_2 &= \frac{x^n_2 + 0x_0 + ax_1 + ax_3}{1 + 2a} \\ x_3 &= \frac{x^n_3 + 0x_0 + 0x_1 + a x_2}{1 + 2a} \end{aligned} $$
We initialize all $x_i = 0$ as our initial guess.
In practice, this is already done at the start of the simulation — at frame 0 we start with all densities at 0, and in subsequent frames, we use the previous frame’s values as our initial guess.
Then we iteratively update each $x_i$ using the most recent values available — this is the Gauss-Seidel idea.
(If we used only old values from the previous iteration, that would be Jacobi.)
A key detail is that we don’t need to sum across the entire row of the matrix. We already know that only the immediate neighbors contribute non-zero terms (this is why Stam enfasises the sparsity nature of the matrix), and we know their coefficient values (−a).
Pseudocode for 1D Diffusion with Gauss-Seidel
for (k = 0; k < maxIterations; k++) {
for (i = 1; i < N; i++) {
x[i] = (x0[i] + a * (x[i-1] + x[i+1])) / (1 + 2*a);
}
}
2D diffusion
The concept generallizes easily to 2D and 3D. The core differences are:
- The discretization of the Laplacian changes to include the 4 (in 2D) or 6 (in 3D) neighboring cells.
- We need to find a way to map our multidimensional indices (i, j) or (i, j, k) into a 1D array so we can fill our coefficient matrix and term vectors.
For the 2D case, we already know the discrete form of the Laplacian in the Backward Euler update rule:
$$ x_0(i,j) = x(i,j) - a \left( x(i-1,j) + x(i+1,j) + x(i,j-1) + x(i,j+1) - 4x(i,j) \right) $$
We can rearrange this to solve for $x(i,j)$:
$$ x(i,j) = \frac{x_0(i,j) + a \left( x(i-1,j) + x(i+1,j) + x(i,j-1) + x(i,j+1) \right)}{1 + 4a} $$
This gives us the update rule used by the Gauss–Seidel method for diffusion, which is exactly what Stam implements in his code:
void diffuse (int N, int b, float *x, float *x0, float diff, float dt) {
int i, j, k;
float a = dt * diff * N * N;
for (k = 0; k < 20; k++) {
for (i = 1; i <= N; i++) {
for (j = 1; j <= N; j++) {
x[IX(i,j)] = (x0[IX(i,j)] + a * (
x[IX(i-1,j)] + x[IX(i+1,j)] +
x[IX(i,j-1)] + x[IX(i,j+1)]
)) / (1 + 4 * a);
}
}
set_bnd(N, b, x);
}
}
That said, I still want to show how the problem looks when expressed in matrix form, so we can see how its sparsity and diagonal dominance extend to multiple dimensions.
Below is the fully expanded $(9\times9)$ matrix system $(A\mathbf{x}=\mathbf{b})$ for a $(3\times3)$ grid (so $(N=n^2=9)$ unknowns). I used the usual row-major indexing:
$$ m = \text{index}(i,j) = i + 3j \quad\text{for}\quad i,j\in{0,1,2} $$
The grid would look something like (each cell is indexed $x_{i,j}^{m}$)
$$ \begin{bmatrix} x_{0,0}^{0} & x_{1,0}^{1} & x_{2,0}^{2} \\ x_{0,1}^{3} & x_{1,1}^{4} & x_{2,1}^{5} \\ x_{0,2}^{6} & x_{1,2}^{7} & x_{2,2}^{8} \end{bmatrix} $$
Matrix A:
$$ \begin{bmatrix} 1+4a & -a & 0 & -a & 0 & 0 & 0 & 0 & 0 \\ -a & 1+4a & -a & 0 & -a & 0 & 0 & 0 & 0 \\ 0 & -a & 1+4a & 0 & 0 & -a & 0 & 0 & 0 \\ -a & 0 & 0 & 1+4a & -a & 0 & -a & 0 & 0 \\ 0 & -a & 0 & -a & 1+4a & -a & 0 & -a & 0 \\ 0 & 0 & -a & 0 & -a & 1+4a & 0 & 0 & -a \\ 0 & 0 & 0 & -a & 0 & 0 & 1+4a & -a & 0 \\ 0 & 0 & 0 & 0 & -a & 0 & -a & 1+4a & -a \\ 0 & 0 & 0 & 0 & 0 & -a & 0 & -a & 1+4a \end{bmatrix} $$
Important note about boundaries:
The matrix shown above assumes that every cell has 4 neighbors, meaning the discrete Laplacian always uses all four surrounding values. Under this assumption, every diagonal entry is $1 + 4a$.
This is exactly how Stam treats the problem in his paper and implementation — the solver itself doesn’t worry about missing neighbors or special boundary handling. Instead, all cells are treated uniformly, and any boundary conditions are enforced afterward through the set_bnd() function, which directly modifies the edge values after each relaxation step.
If we were to remove that assumption — i.e., not include any ghost contributions outside the grid — then boundary cells would have fewer valid neighbors. Corners would have only 2 neighbors, edges 3, and interior cells 4. In that case, each corresponding diagonal term would change to $1 + (\text{number of neighbors})*a$, so corners become $1 + 2a$, edges $1 + 3a$, and only interior cells remain $1 + 4a$.
Closing thoughts
To me, it’s amazing how a few apparently simple lines of code carry so much mathematical depth and careful reasoning.
I’m hoping in the future I get to explore other methods for improving the approximation of the derivative (I’ve heard of Runge-Kutta method which is basically backwards euler but with more steps) and more accurate, faster solvers for the resulting linear systems (a little bird told me that Conjugate Gradient with some tricks is the way to go).
I should also dive deeper into boundary conditions and study how different types (Dirichlet, Neumann, periodic, etc.) influence the simulation.
Another thing Stam doesn’t really touch on in this paper is how many iterations of the Gauss–Seidel relaxation to use. He defaults to 20, and I currently do the same, but ideally, I should measure the residual error at each iteration and determine the minimum number needed for a visually stable and accurate result.
Anyway, those are future adventures. For now, here’s our implicit diffusion solver in action — try bumping up the diffusion coefficient and see how beautifully stable it remains:
Following Velocity
The final term of the density solver is the advection term (sometimes also called convection or transport):
$$ -(\mathbf{u} \cdot \nabla)\rho $$
It describes how the density is transported along the velocity field. In Stam’s equation, for a constant velocity field u and no diffusion or sources, we can write:
$$ \frac{d\rho}{dt} = -(\mathbf{u} \cdot \nabla)\rho $$
or equivalently,
$$ \frac{d\rho}{dt} + (\mathbf{u} \cdot \nabla)\rho = 0 $$
This just means that the density is moving around due to the velocity field, but the total amount of density in the scene stays the same — it’s just being carried along.
Understanding the advection term
The operator $(\mathbf{u} \cdot \nabla)$ takes the gradient of a scalar field and dots it with u.
In multivariable calculus, this is known as the directional derivative — it measures the rate of change of a function along some direction.
While the gradient gives the rate of change in the direction of the coordinate axes, the directional derivative gives the rate of change along an arbitrary vector. In our case, that vector is the velocity field, which literally defines how the density moves in space.
So $(\mathbf{u} \cdot \nabla)\rho$ tells us how the scalar $\rho$ (our density) changes as we move along the flow defined by u. In 2D, it expands as:
$$ (\mathbf{u}\cdot\nabla)\rho = u_x \frac{\partial p}{\partial x} + u_y \frac{\partial \rho}{\partial y}. $$
Deriving an update rule (Forward Euler)
Let’s start simple.
Assume the velocity field u is constant and known (we’re defining it).
We want an update rule that gives us $\rho^{n+1}$ (density at time $t+\Delta t$) from $\rho^n$ (density at time $t$).
Just like in the diffusion step, we can start with the Forward Euler integration scheme:
$$ \rho^{n+1} = \rho^n + \Delta t\, f(\rho^n) $$
where $f(\rho)$ approximates $\frac{d\rho}{dt}$.
Here,
$$ f(\rho) = -(\mathbf{u}\cdot\nabla)\rho. $$
Let’s approximate the spatial derivatives using forward differences (just for now, the simplest possible):
$$ \frac{\partial \rho}{\partial x} \approx \frac{\rho_{i+1,j} - \rho_{i,j}}{h}, \quad \frac{\partial \rho}{\partial y} \approx \frac{\rho_{i,j+1} - \rho_{i,j}}{h}. $$
Plugging into our equation:
$$ \rho_{i,j}^{n+1} = \rho^n_{i,j} - \Delta t \left( u_x \frac{\rho^n_{i+1,j} - \rho^n_{i,j}}{h} + u_y \frac{\rho^n_{i,j+1} - \rho^n_{i,j}}{h} \right). $$
This is a simple, explicit, and easy-to-implement update rule.
However, as we’ve already seen in the diffusion step, explicit schemes can become unstable for large enough time steps.
In this case, the “effective timestep” is related to the CFL number:
$$ \nu = \frac{|\mathbf{u}|\Delta t}{h}. $$
If $\nu > 1$, the update becomes unstable — the simulation starts to blow up or oscillate wildly. (TODO: actually stuy, derive and test the CFL condition for this scheme.)
Trying the backward (implicit) approach
Let’s try something more stable.
We’ll use Backward Euler as the time integration scheme, and a central difference to approximate the spatial derivative.
In 1D, the advection equation becomes:
$$ \frac{\partial \rho}{\partial t} = -u \frac{\partial \rho}{\partial x}. $$
Using Backward Euler in time:
$$ \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = -u \frac{\rho_{i+1}^{n+1} - \rho_{i-1}^{n+1}}{2h}. $$
Rearranging terms:
$$ -\frac{u\Delta t}{2h}\, \rho_{i-1}^{n+1} + \rho^{n+1}i + \frac{u\Delta t}{2h}\, \rho{i+1}^{n+1} = \rho^n_i. $$
Letting $a = \frac{u\Delta t}{2h}$, this becomes:
$$ -a\,\rho^{n+1}_{i-1} + \rho^{n+1}i + a\,\rho^{n+1}{i+1} = \rho^n_i. $$
If we write this for every grid cell, we get a tridiagonal system with off-diagonal terms $(-a, +a)$.
In matrix form (for 5 cells, as an example):
$$ \begin{bmatrix} 1 & a & 0 & 0 & 0\\ -a & 1 & a & 0 & 0\\ 0 & -a & 1 & a & 0\\ 0 & 0 & -a & 1 & a\\ 0 & 0 & 0 & -a & 1 \end{bmatrix} \begin{bmatrix} \rho^{n+1}_1\\ \rho^{n+1}_2\\ \rho^{n+1}_3\\ \rho^{n+1}_4\\ \rho^{n+1}_5 \end{bmatrix} = \begin{bmatrix} \rho^n_1\\ \rho^n_2\\ \rho^n_3\\ \rho^n_4\\ \rho^n_5 \end{bmatrix}. $$
This is solvable, but it’s not ideal.
To ensure Gauss–Seidel convergence, the matrix must either be symmetric positive-definite — which it isn’t, since the off-diagonals have opposite signs and the coefficients depend on velocity (we don’t want to constrain ourselves to constant velocity) — or at least strictly or irreducibly diagonally dominant which implies:
The diagonal entry is $a_{ii} = 1$, and the sum of off-diagonal magnitudes per row is $2|a|$.
For diagonal dominance we need:
$$ |a_{ii}| \ge \sum_{j\neq i} |a_{ij}| $$
which gives:
$$ 1 \ge 2|a| \quad \Rightarrow \quad |u|\Delta t / h \le 0.5 $$
That’s a velocity cap directly tied to the grid resolution and timestep — definitely not something we want to enforce in a fluid sim.
Without it, convergence isn’t guaranteed.
And that’s just in 1D. In 2D or 3D, the situation gets even nastier: the velocity has multiple components ($u_x, u_y, u_z$), so the matrix coefficients vary per cell and direction, breaking both symmetry and uniformity. The matrix stops having a nice constant stencil — every row can look different.
That’s what Stam refers to when he says:
“However, the resulting linear equations would now depend on the velocity, making it trickier to solve.”
He means that since the velocity field changes every frame, the matrix coefficients change too — you’d have to rebuild a new, non-symmetric, potentially non-convergent system at each timestep.
So even though this implicit scheme avoids the usual explicit-time instability, it introduces its own headaches: asymmetric matrices, poor convergence guarantees, and a dependence on a constantly changing velocity field.
Other ideas
Apparently, there are better ways to do this within the Eulerian framework. From the Fluid Simulation for Computer Graphics book:
“In general, biasing a finite difference to the direction that flow is coming from is called upwinding. Most advanced Eulerian schemes are upwind-biased schemes that do this with more accurate finite difference formulas.”
For now, I’m not diving into that. Let’s see how Stam solved this instead — by taking a completely different perspective on the problem.
Before jumping into that, it’s worth clarifying what this means. So far, everything we’ve done lives in the Eulerian viewpoint — we look at fixed points in space and track how quantities like density and velocity change over time at those points. Each grid cell is basically a small “sensor” watching the fluid flow through it.
The Lagrangian viewpoint, on the other hand, treats the fluid as a swarm of individual “particles” that we follow through space and time. Each particle carries its own position, velocity, etc., and we observe how those quantities change as the particle moves.
Semi-Lagrangian Advection
To solve the advection step, Stam takes a page from the Lagrangian viewpoint:
“The key idea behind this new technique is that moving densities would be easy to solve if the density were modeled as a set of particles. In this case we would simply have to trace the particles through the velocity field. For example, we could pretend that each grid cell’s center is a particle and trace it through the velocity field…”
As he states, we can easily write an update rule for our particles’ positions:
$$ \mathbf{p}^{n+1} = \mathbf{p}^n + \Delta t \frac{d\mathbf{p}}{dt} $$
where $\frac{d\mathbf{p}}{dt}$ is the velocity, given by the velocity field itself.
This is why the method is called Semi-Lagrangian — we’re not truly simulating discrete particles; we’re still using a grid, but we think about each cell as if it were a particle moving through the flow.
Tracing Backwards Through the Velocity Field
If we were to trace particles forward, the question would arise of how to turn particles ammounts back into density values. We’d likely have to compute where each one lands and how its mass contributes to nearby cells. Apparently several techniques were developed using that idea, but Stam proposes a simpler alternative: trace backwards.
We know the current density field $\rho^n$. To find the density at a new time step $\rho^{n+1}$, we follow the velocity field backward from each cell center to find where the density came from:
$$ \mathbf{p}^n = \mathbf{p}^{n+1} - \Delta t \, \mathbf{u} $$
This is the key insight: instead of pushing particles forward, we pull density backward from where it came.
for (let i = 1; i <= N; i++) {
for (let j = 1; j <= N; j++) {
let x = i - deltaT * velX[i][j];
let y = j - deltaT * velY[i][j];
// interpolate density at (x, y)
}
}
Note:
This assumes that the grid indices (i, j) correspond to the center of each cell.
If grid coordinates represent cell corners instead, you would need to offset by half a cell width (+h/2) when tracing positions.
Bilinear Interpolation
Once we have the backtraced position $(x, y)$, we convert it to grid indices. The integer part gives us the upper-left cell center (this of course depends on how you construct your grid), while the fractional part tells us the percentage offset we have in relation to nearby cells. We then use bilinear interpolation (weighted average of the four surrounding cells) to sample the previous density field:
(i0, j0) (i0+1, j0)
+--------------+
| |
| (x, y) |
| |
+--------------+
(i0, j0+1) (i0+1, j0+1)
- $i_0 = \lfloor x \rfloor$, $i_1 = i_0 + 1$
- $j_0 = \lfloor y \rfloor$, $j_1 = j_0 + 1$
- $s = x - i_0$, $t = y - j_0$
Interpolate horizontally along the bottom row:
$$ f(x,j_0) = (1 - s)\,\rho_{i_0,j_0} + s\,\rho_{i_1,j_0} $$
Interpolate horizontally along the top row:
$$ f(x,j_1) = (1 - s)\,\rho_{i_0,j_1} + s\,\rho_{i_1,j_1} $$
Interpolate vertically between those two:
$$ \rho(x,y) = (1 - t)\,f(x,j_0) + t\,f(x,j_1) $$
This can be written compactly as:
$$ \rho(x, y) = (1 - s)(1 - t)\rho_{i_0,j_0} + s(1 - t)\rho_{i_1,j_0} + (1 - s)t\rho_{i_0,j_1} + st\rho_{i_1,j_1} $$
On Numerical Diffusion
Here’s where things get tricky. Bilinear interpolation introduces numerical diffusion: sharp features in the density field quickly smooth out, and with coarse grids or large time steps, the density may even dissipate completely.
This diffusion isn’t physical — it’s a numerical artifact of interpolation.
Stam acknowledges this in the paper, and many later works have explored higher-order advection schemes to mitigate the problem, which I hope to explore in the future.
Quick Note on Units and Index Space
We define our domain to have a “physical” size of 1 (you can think of it as 1 meter).
This means $x, y \in [0,1]$.
However, our grid has $N$ cells, so each cell is $1/N$ units wide.
When we write:
$$ x = i - \Delta t_0 \, u_x $$
we are working in index space, not world space.
If our velocity field is expressed in world units (e.g., 0.2 means “move across 20% of the domain per second”), we need to convert it to cells per second by multiplying by $N$.
This is why Stam defines:
$$ \text{dt0} = \Delta t \times N $$
It’s simply converting the velocity from world-space units to cell-space units.
Conclusion
And there we have it — densities moving smoothly along a velocity field!
If you set the diffusion rate to zero, you’ll see how the density purely follows the flow, though over time, interpolation errors cause it to blur and fade.
Evolving Velocities
Now that we can move densities around a velocity field, the next step is to evolve the velocity field itself — to make it change over time according to the physics of fluids.
The behavior of a velocity field is governed by the Navier–Stokes equation:
$$ \frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} $$
Each term describes a different physical process:
External forces - $\mathbf{f}$
This term represents anything that directly pushes or pulls the fluid — gravity, wind, a fan, in our case, the mouse dragging through the simulation grid. It’s the easiest part to understand: you’re simply adding velocity to certain points in the field.
Viscous diffusion - $\nu \nabla^2 \mathbf{u}$
This describes how momentum (velocity) spreads out over time, similar to how heat diffuses through a material.
The constant $\nu$ (nu) is the viscosity coefficient:
- A high viscosity (like honey) means the fluid resists motion — disturbances don’t spread easily.
- A low viscosity (like air or water) means changes in velocity diffuse quickly across the field.
Mathematically, the Laplacian $\nabla^2 \mathbf{u}$ measures how different a velocity is from its neighbors. If a cell’s velocity is much higher than the surrounding ones, this term will act to smooth it out over time.
$$ \frac{\partial \mathbf{u}}{\partial t} = \nu \nabla^2 \mathbf{u} $$
If $\nu$ is large, then $\nu \nabla^2 \mathbf{u}$ is large — meaning the change of velocity over time is strong, and the smoothing/spreading effect of velocity differences happens quickly.
That’s why the velocity field becomes more uniform: gradients and small details disappear faster.
If $\nu$ is small, then $\nu \nabla^2 \mathbf{u}$ is weak — velocity differences persist longer, and the fluid retains sharper local “structures” like vortices.
This corresponds to low-viscosity fluids (like air), where momentum isn’t diffused away as quickly and motion feels more dynamic.
Self-advection — $-(\mathbf{u} \cdot \nabla)\mathbf{u}$
This is the trickiest part conceptually. It says that the velocity field moves itself.
If you imagine releasing a puff of smoke into a wind field, the smoke is carried by the wind — that’s advection. But now imagine the “wind” also carries itself along — faster regions push their momentum into neighboring regions. That’s self-advection.
Physically, this is what gives fluids their swirling, chaotic behavior. Mathematically, it’s what makes the Navier–Stokes equations nonlinear. Read a little more about this in the linear vs nonlinear background chapter.
Recycling Numerical Methods
It’s easy to see that the terms in the velocity equation resemble those in the density equation, so we can reuse most of the same routines we developed earlier to solve it. Stam explicitly states this in his paper and more or less skips re-explaining these parts and goes straight into showing the full velocity step.
void vel_step ( int N, float *u, float *v, float *u0, float *v0,
float visc, float dt )
{
add_source ( N, u, u0, dt );
add_source ( N, v, v0, dt );
SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt );
SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt );
project ( N, u, v, u0, v0 );
SWAP ( u0, u ); SWAP ( v0, v );
advect ( N, 1, u, u0, u0, v0, dt );
advect ( N, 2, v, v0, u0, v0, dt );
project ( N, u, v, u0, v0 );
}
Adding Forces
Adding external forces is handled in the same way as before.
In the density step, add_source() injected more density into a cell.
Here, it injects momentum instead — that is, it modifies the velocity of a cell.
In Stam’s code, this is done via the function get_from_UI(), which writes user input values into u0 and v0.
He doesn’t really explain what happens inside this function, but we can infer that it takes mouse input (like dragging) and translates it into localized velocity additions.
Diffusion and Advection
The next two terms — viscous diffusion and self-advection — use the exact same routines we already wrote for density.
The only distinction is that we’re now dealing with a vector field instead of a scalar field.
By definition, the vector calculus operators (like the Laplacian or directional derivative) act component-wise on vectors.
This means we can compute the x and y components of velocity independently:
-
The Laplacian simply applies to each component:
$$ \nabla^2 \mathbf{u} = (\nabla^2 u_x, \nabla^2 u_y) $$
-
The advection step also updates each component separately, even though both components are used when computing how the velocity field moves.
This is why Stam can just call diffuse() and advect() twice — once for u (the x component) and once for v (the y component).
The structure of the PDE allows this separation because, at this stage, there are no cross-terms that directly mix $u_x$ and $u_y$.
The only new function introduced here is project(), which ensures the resulting velocity field is divergence-free (incompressible) which we’ll have to unpack
Project - Enforcing Imcompressibility
The Incompressibility Condition
When introducing the project() function, Stam writes:
“There is, however, a new routine called
project()which is not present in the density solver. This routine forces the velocity to be mass conserving. This is an important property of real fluids which should be enforced.”
What Stam is referring to here is the fact that, in our simulation, we assume two key physical constraints:
(A) conservation of mass, and
(B) conservation of volume.
From Conservation to Incompressibility
In real life, fluids like air and water can indeed change volume while keeping their mass constant. Those local changes in volume (and therefore density, since $\rho = m/V$) are what give rise to sound waves — tiny oscillations in pressure and density that propagate through the medium.
Fortunately for us, such effects are negligible for visual fluid simulation, so we make the simplifying assumption that density remains constant everywhere.
This means that fluid parcels neither expand nor contract — in other words, they are incompressible.
Under this assumption, the mathematical condition for incompressibility is beautifully simple:
$$ \nabla \cdot \mathbf{u} = 0 $$
This equation states that the divergence of the velocity field is zero everywhere.
Building Intuition
Even without going through the full derivation (which is shown in the background chapter), we can read this equation intuitively.
Imagine a box filled with fluid particles moving along the velocity field.
- If, around some point, the velocity vectors point “more inward” and get shorter as they approach it, that region acts as a sink — it has negative divergence. Particles accumulate there, so the local density increases.
- Conversely, if the velocities point “more outward” and grow in magnitude, the point behaves like a source — it has positive divergence, and the local density decreases as fluid “expands” away.
TODO INCLUDE PICTURES OF THE 3 DIFFERENT CASES OF DIVERGENCE
That behavior describes a compressible fluid — one where regions can locally gain or lose density.
Setting $\nabla \cdot \mathbf{u} = 0$ ensures the opposite: no fluid parcel expands or contracts. The density everywhere stays constant, guaranteeing the flow is incompressible.
Edge cases
It’s important to note that this condition applies at the level of fluid elements — infinitesimal parcels of fluid that move and deform with the flow. Each parcel retains its volume over time, and because these parcels collectively make up the entire domain, enforcing incompressibility locally ensures it globally.
Consider two examples:
- Compressing the box:
Suppose we have fluid inside a sealed box and start shrinking the box. The total amount of fluid (mass) remains constant, but the space it occupies decreases. In an incompressible simulation, this doesn’t cause the fluid to shrink — instead, pressure increases to resist compression, as though the fluid were infinitely stiff. - Injecting more fluid:
If we add more fluid into the box, we’re increasing both the total volume and total mass of the fluid, but as long as the density remains constant, this is still consistent with incompressibility. In practice, we handle this not by modifying the incompressibility equation itself, but by controlling inflow through boundary conditions — specifying where new fluid enters and with what velocity. The interior of the fluid domain still satisfies incompressibility.
Now that we’ve built an intuitive picture of what incompressibility means, let’s look more closely at how it emerges and how the pressure term in project() ensures the condition is satisfied
Deriving the Incompressibility Constraint
Reference: Fluid Simulation for Computer Games, Second Edition - Robert Bridson, Appendix B
This section is taken straight out of Bridson’s work. We begin with the principle of mass conservation. If no fluid is entering or leaving a given region of space there should be no gain/loss of mass along time, the total mass inside that region must remain constant.
Total mass in a region
Let’s consider some fixed region of space, denoted by $\Omega$. So, because $\rho = m / V$ we can think of the total mass in the region as adding up infinitesmly small bits of density multiplied by the tiny volume $d\Omega$, that is integrating the density field $\rho(\mathbf{x}, t)$ over the volume:
$$ M(t) = \iiint_{\Omega} \rho \, d\Omega $$
Rate of change of mass
The rate of change of mass with respect to time only changes as fluid enters or leaves the region. Since mass cannot be created or destroyed inside $\Omega$, the only way $M$ changes is through flux across its boundary $\partial\Omega$. This is basically just Flux from multivariable calculus; which equates to breaking up the surface that bounds $\Omega$ into infinitesimal small pieces of surface $\partial\Omega$ and measuring how much mass is going through a piece which we can do by multiplying density with the dot product between the velocity and surface normal at that “tiny piece”.
Integrating over the whole surface gives:
$$ \frac{dM}{dt} = - \iint_{\partial \Omega} \rho\, (\mathbf{u} \cdot \hat{\mathbf{n}})\, dS $$
Note the integral is negative because $\hat{\mathbf{n}}$ is the outward-pointing normal at the surface, fluid leaving means less mass in region so $M$ decreases (negative rate of change).
Applying the divergence theorem
Then, also from multivariable calculus, we resort to the magnificent divergence theorem which states that adding up all the little bits of outward flow in a volume using a triple integral of divergence gives the total outward flow from that volume, as measured by the flux through its surface.
$$ \iint_{\partial \Omega} \rho\, (\mathbf{u} \cdot \hat{\mathbf{n}})\, dS = \iiint_{\Omega} \nabla \cdot (\rho \mathbf{u})\, d\Omega $$
Substituting into the mass rate equation:
$$ \frac{dM}{dt} = - \iiint_{\Omega} \nabla \cdot (\rho \mathbf{u})\, d\Omega $$
Expressing $dM/dt$ as a volume integral
On the other hand, from the definition of $M(t)$:
$$ M(t) = \iiint_{\Omega} \rho\, d\Omega $$
we can take the time derivative (assuming $\Omega$ is fixed in space) and move the derivative inside the integral:
$$ \frac{dM}{dt} = \iiint_{\Omega} \frac{\partial \rho}{\partial t}, d\Omega $$
(We can safely move the derivative inside because both the region and the limits of integration are constant — only the integrand, the density, changes with time.)
Equating both expressions for $dM/dt$
Now we have two equivalent expressions for the rate of change of mass:
$$ \iiint_{\Omega} \frac{\partial \rho}{\partial t}\, d\Omega = - \iiint_{\Omega} \nabla \cdot (\rho \mathbf{u})\, d\Omega $$
Since this equality must hold for any region $\Omega$, the integrands themselves must be equal.
This gives the continuity equation:
$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 $$
From here we could simply state that for our fluid to be incompressible we want density to remain constant (non-zero and unchanging) by definition. So $\frac{\partial \rho}{\partial t} = 0$ and pull $\rho$ out of $\nabla \cdot (\rho \mathbf{u})$ yielding the incompressible condition $\nabla \cdot \mathbf{u} = 0$
However, this explanation skips over an important idea — how density changes as a fluid moves. To understand that properly, we need to introduce the material derivative.
The Material Derivative
Before moving forward, let’s clarify what we mean when we talk about properties of “moving particles” in a fluid.
In most texts (including Bridson’s), when we refer to a particle of fluid, we’re not talking about a literal molecule, but rather an infinitesimally small “chunk” or “tagged bit” of fluid — small enough to assume uniform properties inside it, but large enough to contain many molecules.
This particle can have a density, temperature, pressure, etc. — those quantities describe the fluid contained in that moving element, not a single molecule.
So when we say “the density of a particle changes,” we mean the density of that small moving volume of fluid — not the density of matter itself (which indeed stays constant for an incompressible fluid).
It might be easier to picture this idea with temperature, suppose you have a pan of water over a fire:
- The fluid could cool down over time as the sun sets (a time-dependent change).
- Or, as the particle moves into a different region of the fluid where the temperature is higher or lower (a spatial change).
Both effects together describe the total rate of change experienced by that moving particle.
That total rate of change is captured mathematically by the material derivative.
Let $f(t, \mathbf{x}(t))$ represent some property of the fluid (like density, temperature, etc.) at time $t$ and position $\mathbf{x}(t) = (x(t), y(t))$, the current position of a moving fluid particle.
By the chain rule:
$$ \frac{Df}{Dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}\frac{dx}{dt} + \frac{\partial f}{\partial y}\frac{dy}{dt} $$
Here $\frac{dx}{dt} = u_x$ and $\frac{dy}{dt} = u_y$ are the components of the velocity vector $\mathbf{u} = (u_x, u_y)$.
We can compactly express this as:
$$ \frac{Df}{Dt} = \frac{\partial f}{\partial t} + \mathbf{u} \cdot \nabla f $$
where $\nabla = (\partial_x, \partial_y)$ is the del (gradient) operator.
This second term $\mathbf{u} \cdot \nabla f$ represents the directional derivative of $f$ along the direction of motion given by velocity $\mathbf{u}$ — how much $f$ changes as the particle moves through space.
In short, the material derivative tells us how a quantity changes for a specific moving bit of fluid, combining both the explicit change in time and the change due to motion.
Finalizing the Incompressibility Constraint
From the principle of mass conservation we arrived at the continuity equation:
$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 $$
We can expand the divergence term using the product rule for divergence:
$$ \nabla \cdot (\rho \mathbf{u}) = \rho\, (\nabla \cdot \mathbf{u}) + \mathbf{u} \cdot \nabla \rho $$
Substituting that into the continuity equation gives:
$$ \frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho + \rho\, (\nabla \cdot \mathbf{u}) = 0 $$
Now notice that the first two terms together are exactly the material derivative of density:
$$ \frac{D\rho}{Dt} = \frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho $$
So we can rewrite the equation as:
$$ \frac{D\rho}{Dt} + \rho\, (\nabla \cdot \mathbf{u}) = 0 $$
Finally, if we assume the fluid is incompressible, that means each fluid parcel maintains constant density as it moves, $\frac{D\rho}{Dt} = 0$, plugging that in:
$$ \rho\, (\nabla \cdot \mathbf{u}) = 0 $$
and since $\rho \neq 0$, we arrive at the incompressibility condition:
$$ \nabla \cdot \mathbf{u} = 0 $$
Enforcing Incompressibility, from Pressure to Projection
_Reference: Projection_method_(fluid*dynamics)*
When Stam writes the equation for velocity he isn’t actually giving the usual Navier-Stokes equations; The incompressible Navier-Stokes equations are usually written as:
Momentum equation:
$$ \frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u}\cdot\nabla)\mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} - \frac{1}{\rho}\nabla p $$
Incompressibility constraint:
$$ \nabla\cdot\mathbf{u} = 0 $$
Compare that to what Stam writes:
$$ \frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u}\cdot\nabla)\mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} $$
He omits both the incompressibility constraint and the term $- \frac{1}{\rho}\nabla p$, although he doesn’t explicitly talk about these in the paper they are in fact solved for and taken into account in the project() function.
In reality, the first three terms of the momentum equation don’t guarantee a divergence free velocity field, the core of the project function is precisely finding the term $- \frac{1}{\rho}\nabla p$ that constraints $u$ to be divergente free. I have yet to derive the Navier-Stokes momentum equation and therefore I don’t really understand where this term comes from and why it looks the way it does. Fortunately for us we can abstract it away and just say we want to find some vector field $w$ that when subtracted from the current velocity $\mathbf{u}^*$ yields a divergent free vector field:
$$ \mathbf{u}^{n+1}=\mathbf{u}^* - w $$
Deriving Projection
To do this, we start of, as Stam does, with something called Helmholtz-Hodge Decomposition following the work written in the book A Mathematical introduction to fluid mechanics (Chorin and Marsden’s) (Got this part from his 1999 Stable Fluids paper not Real-Time Fluid Dynamics for Games). Helmholtz-Hodge Decomposition states that any vector field can be written as the sum of a divergent-free field and the gradient of some scalar valued function:
$$ \mathbf{u}^* = \mathbf{u}^{n+1} + \nabla q $$
Taking the divergence of both sides gives:
$$ \nabla\cdot \mathbf{u}^* = \nabla\cdot(\mathbf{u}^{n+1} + \nabla q) $$
Since divergence distributes over addition and $\nabla\cdot\mathbf{u}^{n+1}=0$ (by definition of incompressibility), we get:
$$ \nabla\cdot \mathbf{u}^* = \nabla\cdot\nabla q = \nabla^2 q $$
This is a Poisson equation for the scalar field $q$.
Once we solve this equation for $q$, we can subtract its gradient from the intermediate velocity:
$$ \mathbf{u}^{n+1} = \mathbf{u}^* - \nabla q $$
By construction, the resulting field $\mathbf{u}^{n+1}$ satisfies $\nabla\cdot\mathbf{u}^{n+1} = 0$.
This process — first computing an intermediate velocity field $\mathbf{u}^*$ ignoring the pressure term, and then projecting it onto the space of divergence-free fields — is known as Chorin’s projection method.
Computing Projection
Now that we’ve seen how pressure enforces incompressibility conceptually, let’s look at how to actually compute it — i.e., how the project() step works in practice.
From Chorin’s projection method, we know we want to find a scalar field $q$ (which corresponds to the pressure term) such that:
$$ \nabla \cdot \mathbf{u}^* = \nabla^2 q $$
Once we solve for $q$, we subtract its gradient from the intermediate velocity field $\mathbf{u}^*$ to get our divergence-free velocity $\mathbf{u}^{n+1}$:
$$ \mathbf{u}^{n+1} = \mathbf{u}^* - \nabla q $$
So the project() function can be divided into two main parts:
- Solve the Poisson equation $\nabla^2 q = \nabla \cdot \mathbf{u}^*$
- Compute $\mathbf{u}^{n+1} = \mathbf{u}^* - \nabla q$
Solving the Poisson Equation
Computing Divergence
Let’s start with the left-hand side. We already have $\mathbf{u}^*$, so we can compute its divergence:
$$ \nabla \cdot \mathbf{u}^{} = \frac{\partial u_{x}^{}}{\partial x} + \frac{\partial u_{y}^{*}}{\partial y} $$
In discrete form, we approximate derivatives using central differences:
$$ \nabla \cdot \mathbf{u}^{} \approx \frac{u^{}(x + h, y) - u^(x - h, y)}{2h} + \frac{v^(x, y + h) - v^(x, y - h)}{2h} $$
Stam’s code:
h = 1.0 / N;
for (i = 1; i <= N; i++) {
for (j = 1; j <= N; j++) {
div[IX(i, j)] = -0.5 * h * (
u[IX(i+1, j)] - u[IX(i-1, j)] +
v[IX(i, j+1)] - v[IX(i, j-1)]
);
}
}
Notice how Stam’s div array doesn’t actually store $\nabla \cdot \mathbf{u}^$ but instead $-h^2\,\nabla \cdot \mathbf{u}^{}$. We can see this because he is multiplying everything by $-0.5$ (same as divinding by $2$ and multiplying by $-1$) and multiplying by $h$ whereas he should be multiplying by $1/h$, and $h = h^2 * 1/h$
He does this so he doesn’t have to go through extra computation in the next steps, as we will see.
Discritizing the laplacian
We’ve already seen how to discretize the Laplacian in the diffusion step:
$$ \nabla^2 q \approx \frac{q(x+h, y) + q(x-h, y) + q(x, y+h) + q(x, y-h) - 4q(x, y)}{h^2} $$
(where $1/h^2 = N*N$)
Given our Poisson equation $\nabla^2 q = \nabla \cdot \mathbf{u}^*$, this gives us a system of linear equations of the form:
$$ Aq = b $$
where:
- $A$ is the matrix representing the discrete Laplacian operator,
- $q$ is the vector of unknowns (the pressure at each grid cell),
- $b$ is the divergence of $\mathbf{u}^*$ we just computed.
Solving the system of linear equations
As previously done when solving the diffusion equation using the Backwards Euler method of integration we can approximate a solution for this linear system of equations by resorting to the Gauss-Seidel iteration method
Solving the discretized Poisson equation for $q(x,y)$ gives us the update rule for the Gauss-Seidel iteration:
Using the discrete Laplacian we’ve seen before, we can write this as:
$$ \frac{ q(x+h, y) + q(x-h, y) + q(x, y+h) + q(x, y-h) - 4q(x, y) }{h^2} = (\nabla \cdot \mathbf{u}^*)(x, y) $$
Let’s rearrange this equation to isolate $q(x, y)$:
$$ -4q(x, y) = h^2 (\nabla \cdot \mathbf{u}^*)(x, y) - (\,q(x+h, y) + q(x-h, y) + q(x, y+h) + q(x, y-h)\,) $$
and then divide both sides by $-4$:
$$ q(x, y) = \frac{ q(x+h, y) + q(x-h, y) + q(x, y+h) + q(x, y-h) - h^2 (\nabla \cdot \mathbf{u}^*)(x, y) }{4} $$
In Stam’s code he just adds div[IX(i, j)] because it’s already storing $- h^2 (\nabla \cdot \mathbf{u}^*)$:
for (k = 0; k < 20; k++) {
for (i = 1; i <= N; i++) {
for (j = 1; j <= N; j++) {
p[IX(i, j)] = (div[IX(i, j)] +
p[IX(i-1, j)] + p[IX(i+1, j)] +
p[IX(i, j-1)] + p[IX(i, j+1)]) / 4;
}
}
set_bnd(N, 0, p);
}
Computing the Divergence-Free Velocity
Finally, we can compute the divergence-free velocity field:
$$ \mathbf{u}^{n+1} = \mathbf{u}^* - \nabla q $$
We again use central differences to approximate the gradient:
$$ \frac{\partial q}{\partial x} \approx \frac{q(x+h, y) - q(x-h, y)}{2h}, \qquad \frac{\partial q}{\partial y} \approx \frac{q(x, y+h) - q(x, y-h)}{2h} $$
Stam implements this step as:
for (i = 1; i <= N; i++) {
for (j = 1; j <= N; j++) {
u[IX(i, j)] -= 0.5 * (p[IX(i+1, j)] - p[IX(i-1, j)]) / h;
v[IX(i, j)] -= 0.5 * (p[IX(i, j+1)] - p[IX(i, j-1)]) / h;
}
}
The result is a velocity field that by construction satisfies the incompressibility condition:
$$ \nabla \cdot \mathbf{u}^{n+1} = 0 $$
And as always, don’t forget to apply boundary conditions after each step.
Projecting Twice??
void vel_step ( int N, float *u, float *v, float *u0, float *v0,
float visc, float dt )
{
add_source ( N, u, u0, dt );
add_source ( N, v, v0, dt );
SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt );
SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt );
project ( N, u, v, u0, v0 );
SWAP ( u0, u ); SWAP ( v0, v );
advect ( N, 1, u, u0, u0, v0, dt );
advect ( N, 2, v, v0, u0, v0, dt );
project ( N, u, v, u0, v0 );
}
Looking at Stam’s code, we can see the projection step is called twice — once after diffusion and once after advection. From our derivation, it’s not immediately obvious why this is necessary, and Stam only briefly says: “We do this because the advect() routine behaves more accurately when the velocity field is mass conserving.”
That’s really the key point: advection assumes the velocity field is divergence-free. When the fluid moves itself around, it needs to do so in a way that doesn’t create or destroy volume. As Bridson puts it, “When we move fluid around and want it to conserve volume, the velocity field we are moving it in must be divergence-free… the sequence of our splitting matters a lot!” So we first project after diffusion to ensure advection starts with a valid, mass-conserving velocity field.
Why not project even earlier, right after the user input? The reason is that it wouldn’t help much — the diffusion step that follows is not self-referenced and would immediately break incompressibility again, wasting computation. So instead, Stam waits until diffusion is done and fixes all accumulated divergence in one go.
Finally, a second projection is needed after advection, since the numerical advection step itself can introduce divergence errors. The result is that we end up projecting exactly twice — the minimal number of times needed to keep the simulation both accurate and efficient. This is my understanding of the subject at the moment, in the future probably would be valuable to plot divergence errors inbetween steps to get some actual data on the subject.
Conclusion
Finally densities moving along a velocity field that’s diffusing and advecting itself and stays (mostly) divergent-free.
The MAC Grid
There are several places where we could start improving our fluid simulator.
I decided to begin with spatial discretization, since it will lead to changes in our data structures and algorithms—and the earlier we make those changes, the better.
Our simulation ultimately solves the incompressible Navier–Stokes equations for velocity and pressure. The density equation is only there to render smoke; the real physics lives in the velocity and pressure solve. So first, a quick recap of the equations we care about (written in the style used in Bridson):
The momentum equation:
$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\cdot \nabla \mathbf{u} + \frac{1}{\rho}\nabla p = \mathbf{g} + \nu\nabla^2 \mathbf{u} $$
Using operator splitting, we solve these subproblems individually:
-
Self-advection:
$$ \frac{\partial \mathbf{u}}{\partial t} = -\mathbf{u}\cdot \nabla \mathbf{u} $$
-
Viscous diffusion:
$$ \frac{\partial \mathbf{u}}{\partial t} = \nu\nabla^2 \mathbf{u} $$
The pressure term comes from the incompressibility condition:
$$ \nabla\cdot\mathbf{u} = 0 $$
This gives a Poisson equation:
$$ \nabla^2 p = \nabla\cdot \mathbf{u} $$
After solving for $p$, we update velocities using $-\nabla p$.
All of this requires accurate numerical approximations of spatial derivatives:
- directional derivatives for advection (though semi-Lagrangian avoids this)
- Laplacians $\nabla^2 u$, $\nabla^2 p$
- divergence $\nabla\cdot u$
- gradient $\nabla p$
We use central differences because they are second-order accurate ($O(h^2)$) and unbiased, unlike forward or backward differences which are only $O(h)$ and directionally one-sided.
But central differences on a collocated grid hide a major problem.
When collocated grids fail
Consider the central difference approximation of the derivative:
$$ \frac{q_{i+1} - q_{i-1}}{2h} $$
Notice something strange: the value at $i$ never appears.
This means certain oscillatory patterns, like $q_i = (-1)^i$, produce a numerical derivative of zero everywhere, even though the function is wildly varying.
Mathematically, the discretization has a non-trivial nullspace: many non-constant fields look like they have zero derivative or zero divergence. When computing divergence in particular, this leads to highly compressible fields being interpreted as perfectly divergence-free, which breaks the pressure projection. The result: fluid volume gain/loss, instability, and general sickness.
One of the cleanest and most widely adopted solutions to this is the MAC grid.
The MAC Grid (Marker-and-Cell)
The idea is simple:
Instead of storing velocities at cell centers, we store them on cell faces.
This gives us a staggered grid:
- Horizontal velocities $u$ live on vertical cell faces.
- Vertical velocities $v$ live on horizontal cell faces.
- Pressure remains at cell centers.
This means we need to store extra rows/columns:
- $u_{i-\frac12,j}$: size $(n_x+1)\times n_y$
- $v_{i,j-\frac12}$: size $n_x\times (n_y+1)$
With this staggering, the central difference for a first derivative aligns perfectly:
$$ \left.\frac{\partial q}{\partial x}\right|i \approx \frac{q{i+\frac12} - q_{i-\frac12}}{h} $$
Note the denominator: $h$ not $2h$. If you expand $q(x\pm h/2)$ in Taylor series, you get this naturally, though it’s easy to see in the image why that is.
Now we must rewrite our update routines to use these staggered velocities—but as we’ll see, most parts barely change.
Diffusion on a staggered grid
Diffusion uses a Laplacian, which remains a 5-point stencil. For the $u$-component:
$$ u^{n+1}{i+\frac12,j} = u^n{i+\frac12,j} + \Delta t,\nu, \frac{ u^{n+1}{i+\frac32,j} + u^{n+1}{i-\frac12,j} + u^{n+1}{i+\frac12,j+1} + u^{n+1}{i+\frac12,j-1} - 4u^{n+1}_{i+\frac12,j} }{h^2} $$
This looks more complicated, but structurally it’s the same diffusion solve—your code for diffusion barely changes.
Projection on a staggered grid
1. Solve the Poisson equation for pressure
The continuous equation is:
$$ \nabla^2 p = \nabla\cdot u $$
Written in discrete form:
$$ \frac{ p_{i-1,j}+p_{i+1,j}+p_{i,j-1}+p_{i,j+1}-4p_{i,j} }{h^2} = \frac{u_{i-\frac12,j} - u_{i+\frac12,j}}{h} + \frac{v_{i,j-\frac12} - v_{i,j+\frac12}}{h} $$
Notice everything lines up perfectly because the velocity components are already sampled at the proper faces.
2. Subtract the pressure gradient to enforce incompressibility
We need $\nabla p$ on faces, but pressure lives at cell centers.
Fortunately, the staggered layout makes this trivial:
$$ \left.\frac{\partial p}{\partial x}\right|{i+\frac12,j} \approx \frac{p{i+1,j} - p_{i,j}}{h} $$
$$ \left.\frac{\partial p}{\partial y}\right|{i,j+\frac12} \approx \frac{p{i,j+1} - p_{i,j}}{h} $$
Then the projection step becomes:
$$ u^{n+1}{i+\frac12,j} = u{i+\frac12,j} - \Delta t, \frac{1}{\rho} \frac{p_{i+1,j} - p_{i,j}}{h} $$
$$ v^{n+1}{i,j+frac12} = v{i,j+\frac12} - \Delta t, \frac{1}{\rho} \frac{p_{i,j+1} - p_{i,j}}{h} $$
That’s all the projection changes: just different indexing.
Semi-Lagrangian advection on a staggered grid
Recall the algorithm on a collocated grid:
- Assume a “particle’’ begins at the cell center.
- Compute its backtraced position using the local velocity.
- Bilinearly interpolate the velocity field at that position.
On a staggered grid, the idea is the same, but we must be more careful about where velocities are defined.
1. Getting a velocity at a cell center
If we’re advecting a cell-centered quantity (like density), we need a velocity at the center.
We compute it by averaging the two velocities around each face:
$$ u_{i,j} = \frac{u_{i-\frac12,j} + u_{i+\frac12,j}}{2} $$
$$ v_{i,j} = \frac{v_{i,j-\frac12} + v_{i,j+\frac12}}{2} $$
Then we backtrace.
2. Advecting face-centered velocities
This is trickier. When advecting $u$-velocities at a face $(i+\frac12,j)$, we need the velocity at that face, not the center.
Bridson gives convenient averaging formulas:
$$ u^{\text{center}}{i,j} = \frac{u{i-\frac12,j} + u_{i+\frac12,j}}{2} $$
$$ v^{\text{center}}{i,j} = \frac{v{i,j-\frac12} + v_{i,j+\frac12}}{2} $$
Face-centered velocities:
$$ u\text{-face at }(i+\tfrac12,j): \quad v = \frac{ v_{i,j-\frac12} + v_{i,j+\frac12} + v_{i+1,j-\frac12} + v_{i+1,j+\frac12} }{4} $$
$$ v\text{-face at }(i,j+\tfrac12): \quad u = \frac{ u_{i-\frac12,j} + u_{i+\frac12,j} + u_{i-\frac12,j+1} + u_{i+\frac12,j+1} }{4} $$
These averages give a consistent, unbiased velocity at the face location.
3. Adjusting the backtrace target
Because velocities live on faces, the “particle’’ representing $u_{i+\frac12,j}$ starts at:
$$ x = (i+\tfrac12)h, \qquad y = jh $$
so after backtracing, we end up at a position that is not aligned to the staggered grid.
Thus we must do:
4. Bilinear interpolation on a staggered grid
Interpolation is the same as before, except the samples come from a grid offset by half a cell.
You bilinearly interpolate on the appropriate staggered array (the $u$-grid or the $v$-grid) using the same formula as usual; only the indexing shifts.
TODO review accuracy with book, insert image and code blocks
Background
If you really had too you could go much further back in math then this. These are just the concepts I really wanted to know in order to understand how the fluid sim stuff works. This guy is amazing, these videos where extremly valuable to making this.
Taylor Series
Suppose we have an unknown function $f(x)$, but we somehow have access to its derivatives at a certain point. We might then try to approximate the function near that point using a polynomial expansion.
The idea is that if our polynomial matches not only the value of $f$ but also some of its derivatives at that point, it can locally mimic the behavior of $f$ quite well.
Polynomial Approximations Around $x = 0$
Let’s start by approximating $f(x)$ around the point $x = 0$.
-
Zeroth-order approximation (degree 0):
$$ p(x) = f(0) $$
This is a constant approximation — clearly very crude, but it matches $f(x)$ at $x = 0$.
-
First-order approximation (degree 1):
We now require that our polynomial $p(x)$ also matches the first derivative of $f$ at 0:
$$ p(0) = f(0), \quad p’(0) = f’(0) $$
The simplest polynomial satisfying these conditions is:
$$ p(x) = f(0) + f’(0)x $$
This gives us the tangent line to $f$ at $x = 0$.
-
Second-order approximation (degree 2):
We now add a quadratic term and require that $p’‘(0) = f’’(0)$:
$$ p(x) = f(0) + f’(0)x + \frac{1}{2}f’’(0)x^2 $$
Checking the derivatives:
$$ p’(x) = f’(0) + f’‘(0)x \Rightarrow p’(0) = f’(0) $$
$$ p’‘(x) = f’‘(0) \Rightarrow p’‘(0) = f’’(0) $$
So this quadratic polynomial matches $f$, $f’$, and $f’’$ at $x = 0$.
The General Pattern
Continuing this process indefinitely gives the Maclaurin series, which is just a Taylor series expanded around $x = 0$:
$$ p(x) = f(0) + f’(0)x + \frac{f’’(0)}{2!}x^2 + \frac{f^{(3)}(0)}{3!}x^3 + \dots + \frac{f^{(n)}(0)}{n!}x^n + \dots $$
This provides an increasingly accurate local approximation of $f(x)$ as more terms are added.
Taylor Series Around an Arbitrary Point $x = c$
If instead we want to approximate $f(x)$ around some point $x = c$, we shift the expansion:
$$ p(x) = f(c) + f’(c)(x - c) + \frac{f’’(c)}{2!}(x - c)^2 + \frac{f^{(3)}(c)}{3!}(x - c)^3 + \dots $$
or more compactly:
$$ p(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(c)}{n!} (x - c)^n $$
This is the Taylor series of $f$ around $x = c$.
It represents how we can reconstruct the local behavior of $f$ using its derivatives, which will be crucial when we approximate derivatives numerically (e.g. in finite-difference schemes).
Numerical Differentiation
Numerical differentiation is a technique for approximating the derivative of a function using finite differences between points, rather than computing the derivative analytically. This approach is essential when an analytical derivative is difficult or impossible to obtain, or when we only have discrete data points. It allows us to estimate the rate of change of a function by approximating the slope of a tangent line using a small step size $h$.
The most common finite difference methods are:
- Forward difference: uses the function value at a point and a point ahead.
- Backward difference: uses the function value at a point and a point behind.
- Central difference: averages forward and backward differences for higher accuracy.
These methods differ in their accuracy and how the error scales with the step size $h$.
Approximating the Derivative
Suppose we have a function $f(x)$ for which we cannot compute the derivative analytically, or we are working with a discrete set of data points. How can we approximate $f’(x)$?
Starting from the definition of the derivative:
$$ f’(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}, $$
we can replace the limit with a small, finite $h$ to obtain an approximation. This leads to the forward difference formula:
$$ f’(x) \approx \frac{f(x+h)-f(x)}{h}. $$
Similarly, we can look backward instead of forward:
$$ f’(x) \approx \frac{f(x) - f(x-h)}{h}, $$
which gives the backward difference formula.
We can also combine both forward and backward differences to improve accuracy, giving the central difference: Now, looking at these two, one might think “Why not combine them? Maybe if we average the forward and backward slopes, we get something better.” And indeed, that leads us to the central difference:
$$ f’(x) \approx \frac{f(x+h) - f(x-h)}{2h}. $$
I guess one could come up with a bunch of ways to approximate the derivative at a point, but maybe more importantly we should be thinking how good is any given approximation?
TODO GRAPH THESE 3 BAD BOYS
Understanding the Error
To analyze how accurate these approximations are, we can use a Taylor series expansion.
The Taylor expansion of $f(x+h)$ around $x$ is:
$$ f(x+h) = f(x) + h f’(x) + \frac{h^2}{2} f’’(x) + \frac{h^3}{6} f^{(3)}(x) + \cdots $$
where $f^{(n)}(x)$ is the $n$-th derivative of $f$ at $x$, and the terms continue indefinitely. The notation $o(h^n)$ refers to terms that scale like $h^n$ or smaller, which become negligible as $h \to 0$.
Forward Difference Error
Substitute the Taylor expansion into the forward difference formula:
$$ \frac{f(x+h) - f(x)}{h}. $$
$$ \frac{f(x+h) - f(x)}{h} = \color{aquamarine}{f’(x)} \color{#ff7faa}+ {\frac{h}{2} f’’(x) + \frac{h^2}{6} f^{(3)}(x) + \cdots} $$
- The first term $\color{aquamarine}{f’(x)}$ is exactly what we want: the derivative at $x$.
- The remaining terms $\color{#ff7faa}{\frac{h}{2} f’’(x) + \frac{h^2}{6} f^{(3)}(x) + \cdots}$ are the error introduced by the finite difference approximation.
We usually summarize this error using the leading order term, which is $\mathcal{O}(h)$ here. Intuitively, this means the discrepancy between the true derivative and our approximation scales with $h$.
In other words:
- If you want the error to be 10 times smaller, you roughly need to make $h$ 10 times smaller.
- The higher-order terms (like $h^2, h^3, \dots$) quickly become negligible as $h$ gets small.
- Essentially, $h$ “dominates” the error. For example: if $h = 0.1$, then $h^2 = 0.01$ and $h^3 = 0.001$, so the first term in the error really dictates how accurate our approximation is.
Backward Difference Error
Similarly, for the backward difference:
$$ f’(x) \approx \frac{f(x) - f(x-h)}{h}, $$
expand $f(x-h)$ using Taylor:
$$ f(x-h) = f(x) - h f’(x) + \frac{h^2}{2} f’’(x) - \frac{h^3}{6} f^{(3)}(x) + \cdots $$
Then:
$$ \frac{f(x) - f(x-h)}{h} = \frac{f(x) - \big(f(x) - h f’(x) + \frac{h^2}{2} f’‘(x) - \cdots \big)}{h} = \color{aquamarine}f’(x) \color{#ff7faa}- \frac{h}{2} f’’(x) + \frac{h^2}{6} f^{(3)}(x) - \cdots $$
- The leading order error is again $\mathcal{O}(h)$, but with opposite sign compared to the forward difference.
A Better Way to Think About the Central Difference
So far, we saw that both the forward and backward difference approximations introduce an error term that comes from the second derivative in the Taylor expansion.
That might get you wondering — is there a way to make that error smaller?
Well, since we know where the error comes from (the second derivative term), maybe we can arrange our approximations in such a way that those error terms cancel each other out.
Let’s look again at the Taylor expansions of $f(x+h)$ and $f(x-h)$:
$$ \begin{aligned} f(x+h) &= f(x) + h f’(x) + \frac{h^2}{2!} f’’(x) + \frac{h^3}{3!} f^{(3)}(x) + \cdots,\end{aligned} $$
$$ \begin{aligned} f(x-h) &= f(x) - h f’(x) + \frac{h^2}{2!} f’’(x) - \frac{h^3}{3!} f^{(3)}(x) + \cdots. \end{aligned} $$
Notice how the even-order terms (those with $f’‘(x), f^{(4)}(x), \dots$) have the same sign, while the odd-order terms (those with $f’(x), f^{(3)}(x), \dots$) flip signs.
Now, if we subtract these two equations, the even terms (and thus our second derivative error) cancel out perfectly:
$$ f(x+h) - f(x-h) = 2h f’(x) + \frac{2h^3}{3!} f^{(3)}(x) + \cdots $$
Dividing both sides by $2h$ gives:
$$ \frac{f(x+h) - f(x-h)}{2h} = \color{aquamarine}f’(x) \color{#ff7faa}+ \frac{h^2}{3!} f^{(3)}(x) + \cdots $$
There it is — the central difference formula.
$$ \boxed{f’(x) \approx \frac{f(x+h) - f(x-h)}{2h}} $$
The key insight here is that we didn’t just “average” two estimates.
We canceled their dominant errors by exploiting symmetry.
This gives us a much more accurate approximation — the leading error term is now $\mathcal{O}(h^2)$.
That means if you make $h$ 10× smaller, your error gets roughly 100× smaller!
This same logic — using Taylor expansions to cancel higher-order terms - is the foundation of more advanced finite difference schemes, if you’d like to explore that, look up higher-order finite difference formulas
Numerical Differentiation of the Second Derivative
What about the second derivative?
Formally, we can define it as:
$$ f’‘(x) = \frac{f’(x+h) - f’(x)}{h}. $$
At first glance, you might think: “Okay, I already have formulas for the first derivative — I’ll just plug those in here.” That would probably work, but it tends to introduce more numerical error? (though honestly i haven’t checked this, it might be a fun exerciset)
Instead, let’s stick with the idea we used before — start from the Taylor series and see if we can combine terms to isolate the the second derivative while canceling as much of the rest as possible.
Looking back at the Taylor expansions for $f(x + h)$ and $f(x - h)$:
$$ \begin{aligned} f(x+h) &= f(x) + h f’(x) + \frac{h^2}{2!} f’’(x) + \frac{h^3}{3!} f^{(3)}(x) + \cdots,\end{aligned} $$
$$ \begin{aligned} f(x-h) &= f(x) - h f’(x) + \frac{h^2}{2!} f’’(x) - \frac{h^3}{3!} f^{(3)}(x) + \cdots. \end{aligned} $$
We can see that if we add these two equations together the first derivative and all other odd terms cancel out, leaving only even-order terms:
$$ f(x + h) + f(x - h) = 2f(x) + h^2 f’’(x) + \frac{h^4}{12} f^{(4)}(x) + \mathcal{O}(h^6). $$
So to finish isolating the second derivative we just have to subtract $2f(x)$ and divide by $h^2$
$$ f’’(x) = \frac{f(x + h) - 2f(x) + f(x - h)}{h^2} + \mathcal{O}(h^2). $$
And there we have it — the central difference formula for the second derivative.
Numerical roundoff error
Given all of this, one might think: “Well, we can just make $h$ really, really small and get a super accurate derivative!”
But one would be wrong — because computers can’t store all decimal numbers exactly.
In double-precision floating point format (IEEE 754), we use 64 bits of memory, divided as follows:
- 1 bit for the sign,
- 11 bits for the exponent, and
- 52 bits for the mantissa (plus one implicit bit).
This gives us about $53 \times \log_{10}(2) \approx 15.955$ digits of precision — roughly 16 decimal digits.
That means that for a number with an infinite fractional part, such as $\sqrt{2}$, adding anything smaller than $10^{-16}$ won’t change its stored value:
$$ \sqrt{2} + 10^{-16}/2 = \sqrt{2}. $$
This is important because when performing millions of operations, these tiny rounding errors accumulate.
Let’s look at the central difference formula again and account for rounding error.
Each sample of $f$ will introduce a small roundoff error $\varepsilon_r$.
$$ \frac{df}{dx} \approx \frac{f(x+h) - f(x-h) + 2\varepsilon_r}{2h} + \mathcal{O}(h^2) $$
If we separate the error contributions, the total error $E$ can be approximated as
$$ E \leq \frac{\varepsilon_r}{h} + \frac{m}{6}h^2, $$
where:
- $\varepsilon_r$ is the roundoff error,
- the $\frac{m}{6}h^2$ term comes from the Taylor expansion ($3! = 6$).
- $m = \max |f^{(3)}(x)|$ for $x \in [x-h, x+h]$ (a constant bounding the third derivative),
Notice what happens as $h$ changes:
- As $h \to 0$, the roundoff term $\varepsilon_r / h$ → ∞.
- As $h \to \infty$, the truncation term $(m/6)h^2$ → ∞.
So there must be a sweet spot — a minimum total error at some optimal $h$.
To find it, we minimize $E(h)$:
$$ \frac{dE}{dh} = 0 \quad \Rightarrow \quad -\frac{\varepsilon_r}{h^2} + \frac{m}{3}h = 0 $$
$$ \Rightarrow \quad h_{\text{opt}} = \sqrt[3]{\frac{3\varepsilon_r}{m}}. $$
If we assume $m$ is roughly of order 1 (not too large or small), and $\varepsilon_r \approx 10^{-16}$ for double precision, then:
$$ h_{\text{opt}} \approx \sqrt[3]{3 \times 10^{-16}} \approx 10^{-5}. $$
So in practice, an $h$ value around 10⁻⁵ is near-optimal for central differences in double precision.
Of course, in more advanced or sensitive simulations, computing this value adaptively can be worthwhile.
Laplacian Operator
Numerical Integration of ODEs
Numerical Solvers for Linear Systems
Reference: GeeksForGeeks – Jacobi Method
Given a system of linear equations, i.e.:
$$ \begin{cases} a_{11}x_1 + a_{12}x_2 = b_1 \\ a_{21}x_1 + a_{22}x_2 = b_2 \end{cases} $$
or in matrix form:
$$ A\mathbf{x} = \mathbf{b} \quad \text{where} \quad A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix} $$
We know how to solve these directly using methods like Gaussian elimination. Unfortunately, those methods tend to be computationally expensive for large systems, so very smart people came up with iterative methods that approximate the solution efficiently.
Jacobi Iterative Method
A specific implementation of the Jacobi method requires the following assumptions:
-
The system of linear equations can be written in the form $A\mathbf{x} = \mathbf{b}$, where:
- $A$ is the coefficients matrix
- $\mathbf{x}$ is the vector of unknowns
- $\mathbf{b}$ is the vector of constants
-
The system of linear equations has a unique solution.
Why? Because if multiple solutions exist, there’s no single point for the method to “converge” to — it could just oscillate between possibilities. -
The coefficient matrix $A$ has no zeroes on its main diagonal.
Why? Because each iteration involves dividing by the diagonal element $a_{ii}$. If any of these are zero, the formula breaks. -
The system is diagonally dominant, meaning:
$$ |a_{ii}| \ge \sum_{j \ne i} |a_{ij}| $$
for every row, and for at least one row the inequality is strict.
Why? This ensures the method converges — roughly speaking, it means each variable “mostly depends on itself” rather than its neighbors, which prevents the iterations from diverging or oscillating.
The Algorithm
-
Rewrite each equation to isolate its variable:
$$ x_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j \ne i} a_{ij}x_j \right) $$
-
Make an initial guess for $\mathbf{x}$, typically using zeroes unless we have better prior information.
-
Compute the first approximation $\mathbf{x}^{(1)}$ by substituting the initial guess into the rewritten equations.
-
Repeat step 3 using the most recent approximation to get the next one:
$$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left(b_i - \sum_{j \ne i} a_{ij}x_j^{(k)} \right) $$
until the error between iterations becomes small enough.
Example
Consider the system of linear equations:
$$ \begin{cases} 2x + y + z = 6 \\ x + 3y - z = 0 \\ -x + y + 2z = 3 \end{cases} $$
Checking if the equations are diagonally dominant:
$$ |2| \ge |1| + |1| = 2 \quad \text{(true)} $$
$$ |3| \ge |1| + |-1| = 2 \quad \text{(true)} $$
$$ |2| \ge |-1| + |1| = 2 \quad \text{(true)} $$
Rewriting the system in Jacobi form:
$$ \begin{cases} x = \frac{6 - y - z}{2} \\ y = \frac{-x + z}{3} \\ z = \frac{3 + x - y}{2} \end{cases} $$
Starting with the initial guess $\mathbf{p}^{(0)} = (0, 0, 0)$:
Iteration 1
$$ x^{(1)} = \frac{6 - 0 - 0}{2} = 3 $$
$$ y^{(1)} = \frac{-0 + 0}{3} = 0 $$
$$ z^{(1)} = \frac{3 + 0 - 0}{2} = 1.5 $$
Iteration 2, using $\mathbf{p}^{(1)} = (3, 0, 1.5)$:
$$ x^{(2)} = \frac{6 - 0 - 1.5}{2} = 2.25 $$
$$ y^{(2)} = \frac{-3 + 1.5}{3} = -0.5 $$
$$ z^{(2)} = \frac{3 + 3 - 0}{2} = 3 $$
Iteration 3, using $\mathbf{p}^{(2)} = (2.25, -0.5, 3)$:
$$ x^{(3)} = \frac{6 + 0.5 - 3}{2} = 1.75 $$
$$ y^{(3)} = \frac{-2.25 + 3}{3} = 0.25 $$
$$ z^{(3)} = \frac{3 + 2.25 + 0.5}{2} = 2.875 $$
Iteration 4, using $\mathbf{p}^{(3)} = (1.75, 0.25, 2.875)$:
$$ x^{(4)} = \frac{6 - 0.25 - 2.875}{2} = 1.4375 $$
$$ y^{(4)} = \frac{0 - 1.75 + 2.875}{3} = 0.375 $$
$$ z^{(4)} = \frac{3 + 1.75 - 0.25}{2} = 2.25 $$
After more iterations, the results will get closer to the exact solution.
We can check the error by plugging the results into the equations and seeing how close the residuals are to 0. At iteration 4:
$$ x = \frac{6 - y - z}{2} \Rightarrow 0 = \frac{6 - y - z}{2} - x \Rightarrow \frac{6 - 0.375 - 2.25}{2} - 1.4375 = 0.25 $$
$$ y = \frac{-x + z}{3} \Rightarrow 0 = \frac{-x + z}{3} - y \Rightarrow \frac{-1.4375 + 2.25}{3} - 0.375 = -0.0625 $$
$$ z = \frac{3 + x - y}{2} \Rightarrow 0 = \frac{3 + x - y}{2} - z \Rightarrow \frac{3 + 1.4375 - 0.375}{2} - 2.25 = 0.03125 $$
Jacobi Iteration Visualizer
Gauss-Seidel Iterative Method
The Gauss–Seidel method builds upon the Jacobi method.
The core idea is simple: within each iteration, we immediately use any newly computed values as soon as they are available, instead of waiting for the next iteration.
In the Jacobi method, all updates are done based on the previous iteration $x^{(k)}$:
$$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \ne i} a_{ij} x_j^{(k)} \right) $$
In the Gauss–Seidel method, however, as soon as we compute a new $x_i^{(k+1)}$, we reuse it in the same iteration for computing later components $x_j^{(k+1)}$ (where $j > i$):
$$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \right) $$
This usually yields faster convergence because each iteration immediately benefits from the most recent updates.
Example
Consider again the system:
$$ \begin{cases} 2x + y + z = 6 \\ x + 3y - z = 0 \\ -x + y + 2z = 3 \end{cases} $$
Rewriting in Gauss–Seidel form:
$$ \begin{aligned} x &= \frac{6 - y - z}{2} \ y &= \frac{0 - x + z}{3} \ z &= \frac{3 + x - y}{2} \end{aligned} $$
Starting with the initial guess $(x, y, z) = (0, 0, 0)$:
| Iteration | $x^{(k+1)}$ | $y^{(k+1)}$ | $z^{(k+1)}$ |
|---|---|---|---|
| 1 | $x = 3.000$ | $y = (-3 + 0)/3 = -1.000$ | $z = (3 + 3 - (-1))/2 = 3.5$ |
| 2 | $x = (6 - (-1) - 3.5)/2 = 1.75$ | $y = (-1.75 + 3.5)/3 = 0.583$ | $z = (3 + 1.75 - 0.583)/2 = 2.083$ |
| 3 | $x = (6 - 0.583 - 2.083)/2 = 1.667$ | $y = (-1.667 + 2.083)/3 = 0.139$ | $z = (3 + 1.667 - 0.139)/2 = 2.264$ |
| 4 | $x = (6 - 0.139 - 2.264)/2 = 1.798$ | $y = (-1.798 + 2.264)/3 = 0.155$ | $z = (3 + 1.798 - 0.155)/2 = 2.321$ |
Gauss–Seidel Iteration Visualizer
Linear vs. Nonlinear Equations
The key feature that makes a PDE nonlinear is self-reference — the dependent variable appears multiplied or composed with itself.
If we unpack the advection equation
$$ \frac{\partial u}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u}, $$
we see that the rate of change of $\mathbf{u}$ with respect to time depends on $\mathbf{u}$ itself and its gradient.
This creates a feedback loop: the velocity field modifies itself over time.
In contrast, consider the diffusion term
$$ \frac{\partial u}{\partial t} = \nu \nabla^2 u. $$
If we scale $u$ by a constant factor $a$, the entire equation still holds with the same scaling:
$$ L(a u) = a L(u), $$
which means the relationship is linear.
Superposition Principle
For linear differential equations, superposition holds.
If $u_1$ and $u_2$ are both solutions, then any linear combination
$$ \alpha u_1 + \beta u_2 $$
is also a valid solution.
Linear operators — such as gradient ($\nabla$), divergence ($\nabla \cdot$), Laplacian ($\nabla^2$), or constant multiplication — all satisfy this property:
$$ L(au_1 + bu_2) = aL(u_1) + bL(u_2) $$
For example:
$$ L(u) = \frac{\partial u}{\partial x} \Rightarrow L(au_1 + bu_2) = a\frac{\partial u_1}{\partial x} + b\frac{\partial u_2}{\partial x}. $$
Nonlinear Operator Example
Now let’s take a nonlinear operator:
$$ L(u) = u \frac{du}{dx}. $$
If this were linear, we should have
$$ L(au_1 + bu_2) = aL(u_1) + bL(u_2), $$
but let’s check:
$$ L(au_1 + bu_2) = (au_1 + bu_2) \frac{d(au_1 + bu_2)}{dx} = (au_1 + bu_2)(a\frac{du_1}{dx} + b\frac{du_2}{dx}). $$
Expanding:
$$ a^2u_1\frac{du_1}{dx} + ab\,u_1\frac{du_2}{dx} + ab\,u_2\frac{du_1}{dx} + b^2u_2\frac{du_2}{dx}. $$
Compare that to:
$$ aL(u_1) + bL(u_2) = a,u_1\frac{du_1}{dx} + b,u_2\frac{du_2}{dx}. $$
They are not equal — extra cross terms appear.
This violates superposition, confirming that $L(u) = u \frac{du}{dx}$ is a nonlinear operator.
TODO explain further stuff