Motivation

In the spring of 2026 I started to work on data-driven discovery of PDEs as part of my PhD project, with a focus on equivariant methods for PDE discovery on flat planes, spheres, and hyperplanes. Instead of deriving equations from physical principles the idea is to start from data and try to deduce the governing equations. This is fairly well studied, but the question of how to handle the same problem when symmetries are known (or suspected) or when the domain has non-zero curvature was still unexplored.

Erlend (my doctoral advisor) and I had already worked on the foundations of this problem from the theoretical side, by studying a complete classification of equivariant polynomial differential operators ¹, and we wanted to now apply our results to a practical problem. I started by designing a research plan which involved picking some famous PDEs to test, implementing a numerical integrator to simulate them, and applying our methods against known benchmarks. Doing a literature review, one of the PDEs I encountered was the Kuramoto–Sivashinsky (KS) equation ²,³, and I was immediately fascinated by it!

The KS equation is one of the simplest partial differential equations that exhibits several regimes among which some are spatiotemporal chaos. It is a canonical testbed for studying chaos in extended systems, and it appears in surprisingly many physical contexts.

What struck me most was the dramatic change in behaviour (the different patterns) as a single parameter (the domain size \(L\)) is varied. For small domains the equation is boring as everything decays to zero, but for slightly larger domains, steady cellular patterns emerge (hexagons, rolls, or oscillating cells). Push \(L\) further, and the cells break apart into fully developed spatiotemporal chaos. Three qualitatively different regimes from one parameter. I wanted to see them, and I wanted to play with them. I then needed to write my own simulator to do that, and I decided to build it as a browser-based interactive tool for everyone to enjoy.

This post walks through the equation, the three regimes, and briefly how the simulator works under the hood. You can play with it yourself here.

---

The equation in 1D

Before diving into two dimensions, it's worth looking at the one-dimensional version. In 1D the KS equation is usually written as

\[\partial_t u + \partial_x^2 u + \partial_x^4 u + \frac{1}{2}(\partial_x u)^2 = 0,\qquad x \in [0, L], \quad \text{periodic}.\]

To understand what each term does, we switch to Fourier space (where \(\partial_x \to i k\), turning derivatives into multiplication). Expanding \(u(x,t) = \sum_k \hat{u}_k(t) e^{i k x}\) with \(k = 2\pi j / L\) for \(j \in \mathbb{Z}\), the equation becomes

\[\partial_t \hat{u}_k = (k^2 - k^4)\hat{u}_k + \hat{N}(u)_k,\]

where the linear part breaks down as follows:

Physical space	Fourier space	Effect
\(-\partial_x^2 u\)	\(k^2 \hat{u}_k\)	anti-diffusion (amplifies)
\(-\partial_x^4 u\)	\(-k^4 \hat{u}_k\)	hyperdiffusion (damps)

The two linear terms compete: \(k^2\) wants to grow the mode, \(-k^4\) wants to kill it. Modes with \(\vert k\vert < 1\) grow (anti-diffusion wins), modes with \(\vert k\vert > 1\) decay (hyperdiffusion wins). The nonlinear term \(\hat{N}(u)_k = \mathcal{F}[-\tfrac{1}{2}(\partial_x u)^2]\) couples different \(k\) together, transferring energy between scales.

The classic way to visualize the 1D KS equation is through a spacetime plot: a horizontal slice of the solution stacked vertically over time.

---

The equation in 2D

In two spatial dimensions the KS equation takes the form

\[\partial_t u = -\nabla^2 u - \nu \nabla^4 u - \frac{1}{2}\vert\nabla u\vert^2,\qquad(x,y) \in [0, L]^2, \quad \text{periodic}.\]

The parameter \(\nu\) controls the strength of the hyperdiffusion (in the 1D version above, \(\nu\) was implicitly set to 1).

To analyze the equation, we expand the solution in a 2D Fourier series. Because the domain is periodic, we can write

\[u(x,y,t) = \sum_{\mathbf{k}} \hat{u}_{\mathbf{k}}(t)\, e^{i(k_x x + k_y y)},\]

where the index of the Fourier coefficients \(\mathbf{k} = (k_x, k_y)\) is called wavevector, which label each Fourier mode. The magnitude \(k = \vert\mathbf{k}\vert = \sqrt{k_x^2 + k_y^2}\) is called the wavenumber. A mode with large \(k\) varies rapidly in space (small ripples), while a mode with small \(k\) varies slowly (broad swells). The coefficient \(\hat{u}_{\mathbf{k}}(t)\) is a complex number encoding the amplitude and phase of the mode with wavevector \(\mathbf{k}\) at time \(t\).

Working in Fourier space turns derivatives into multiplication, and just as \(\partial_x\) turns into \(i k_x\) in 1D, in 2D the gradient becomes \(i\mathbf{k} = (i k_x, i k_y)\), and the Laplacian becomes \(-\vert\mathbf{k}\vert^2 = -(k_x^2 + k_y^2)\). Applying this to the KS equation, we split it into a linear part \((-\Delta u - \nu\Delta^2 u)\) and a nonlinear part \((-\tfrac{1}{2}\vert\nabla u\vert^2)\), and transform each:

• The linear part becomes \(\bigl(\vert\mathbf{k}\vert^2 - \nu\vert\mathbf{k}\vert^4\bigr)\,\hat{u}_{\mathbf{k}}\). Crucially, this is diagonal: each Fourier mode \(\hat{u}_{\mathbf{k}}\) is multiplied only by a scalar that depends on \(\mathbf{k}\) alone. There is no mixing between different \(\mathbf{k}\) — mode \(\mathbf{k}\) doesn't talk to mode \(\mathbf{k}'\) through the linear terms.
• The nonlinear part, which we write as \(\hat{N}(u)_{\mathbf{k}} = \mathcal{F}\bigl[-\tfrac{1}{2}\vert\nabla u\vert^2\bigr]\), is not diagonal: it couples different Fourier modes through convolution. This coupling is what transfers energy between scales and ultimately generates chaos.

So in Fourier space the PDE becomes a system of ordinary differential equations, one for each \(\mathbf{k}\):

\[\partial_t \hat{u}_{\mathbf{k}} = \bigl(\vert\mathbf{k}\vert^2 - \nu\vert\mathbf{k}\vert^4\bigr)\,\hat{u}_{\mathbf{k}} \;+\; \hat{N}(u)_{\mathbf{k}}.\]

The linear coefficient depends only on \(k = \vert\mathbf{k}\vert\), so the linear part is diagonal — each mode evolves independently. Because it is diagonal, the linear part can be solved exactly: over a time step \(h\), each mode simply multiplies by

\[\hat{u}_{\mathbf{k}}(t + h) = e^{(k^2 - \nu k^4)h} \,\hat{u}_{\mathbf{k}}(t), \qquad\text{where } k = \vert\mathbf{k}\vert.\]

Just as in 1D, the growth rate \(k^2 - \nu k^4\) pits anti-diffusion (\(k^2\), wants to grow) against hyperdiffusion (\(-\nu k^4\), wants to decay). The two balance at the critical wavenumber

\[k_{\text{crit}} = \frac{1}{\sqrt{\nu}}.\]

Modes with \(k < k_{\text{crit}}\) grow; modes with \(k > k_{\text{crit}}\) decay.

What does this have to do with the domain size? On a periodic domain of length \(L\), the smallest nonzero wavenumber is \(k_{\min} = 2\pi/L\) and you can't fit a wave longer than the box in the box itself. If even this longest wave decays (\(k_{\min} > k_{\text{crit}}\)), the solution dies out. The threshold \(k_{\min} = k_{\text{crit}}\) gives the critical domain size

\[L^\ast = 2\pi\sqrt{\nu}.\]

---

Three regimes, one equation

The ratio \(L / L^\ast\) is the single parameter that determines the qualitative behaviour of the system. Here's what happens as you increase it.

Stable (\(L < L^\ast\))

When the domain is smaller than the critical length, no Fourier mode fits within the unstable band. Every mode decays exponentially, and any initial condition eventually flattens to \(u \equiv 0\). Not terribly exciting, but a good sanity check for the simulator.

Cellular (\(L^\ast < L \lesssim 3L^\ast\))

A small set of modes — those whose wavevectors fit inside the unstable annulus \(0 < \vert\mathbf{k}\vert < 1/\sqrt{\nu}\) — begin to grow. Nonlinear saturation locks them into steady or periodically oscillating patterns: rolls, hexagons, or mixed states depending on the initial condition and the precise value of \(L\).

Chaotic (\(L \gg L^\ast\))

When the domain supports many unstable modes, the nonlinear term transfers energy between them continuously and the patterns never settle. They twist, merge, split, and reform in an endless dance which is spatiotemporal chaos. The field at any given point fluctuates unpredictably and produces a result that is visually mesmerising.

---

How the simulator works

The simulator runs entirely in the browser with the use of a JavaScript snippet. The numerical engine is implemented from scratch following standard methods from the literature, with some tweaks to balance accuracy and performance.

Halfway through the project I came across Ricky Reusser's beautiful WebGL notebook which implements a similar simulator, which I recommend checking out. Especially his WebGPU implementation which is incredibly fast.

Here's a high-level overview of the key components of the simulator I built. I am not an expert in numerical PDEs, so a lot of the design choices were made through trial and error and relying on papers I scavenged from the internet. If you have suggestions for improvements, please let me know!

Spectral discretization

The domain \([0, L]^2\) is discretised into an \(M \times M\) grid. Since the boundary conditions are periodic, the natural basis is the Fourier basis. The solution is represented by its Fourier coefficients \(\hat{u}_{\mathbf{k}}\) on a grid of wavenumbers

\[\mathbf{k} \in \frac{2\pi}{L}\big\{(j_x, j_y) \;\big\vert\; j_x, j_y = 0, \pm 1, \dots, \pm M/2\big\}.\]

Spatial derivatives are computed in Fourier space (multiplication by \(i k_x\) or \(i k_y\)), then transformed back to physical space when needed. The 2D FFT is the workhorse — specifically a hand-written radix-2 Cooley–Tukey implementation. It's not as fast as FFTW, but it runs in a Web Worker so the UI stays responsive.

The nonlinear term

The nonlinear term \(\hat{N}(u) = \mathcal{F}[-\frac{1}{2}\vert\nabla u\vert^2]\) is computed pseudo-spectrally:

1. Compute the physical-space gradient \(\nabla u = \mathcal{F}^{-1}[i\mathbf{k}\,\hat{u}]\)
2. Form \(-\frac{1}{2}\vert\nabla u\vert^2\) pointwise in physical space
3. Transform back: \(\hat{N} = \mathcal{F}[-\frac{1}{2}\vert\nabla u\vert^2]\)

Apparently, there can be an issue of aliasing, which can be addressed by the 2/3 de-aliasing rule: after each step, zero out all Fourier modes with \(\vert j_x\vert > M/3\) or \(\vert j_y\vert > M/3\). This removes the aliasing error at the cost of effectively running on a coarser grid (and about a 2× slowdown). In practice, I don't really see much difference.

Exponential Time Differencing (ETD)

The hyperdiffusion term \(-\nu\nabla^4 u\) is extremely stiff. A stiff term is one that involves drastically different time scales within the same problem and here, the high-wavenumber modes decay orders of magnitude faster than the low-wavenumber ones. Any explicit time-stepping method must take steps small enough to resolve the fastest scale, even if we only care about the slow ones.

Concretely, let's see what happens if we try to use forward Euler. For the highest resolved wavenumber \(\vert\mathbf{k}\vert_{\max} \approx \pi M/L\) the hyperdiffusion term \(-\nu\vert\mathbf{k}\vert^4 \hat{u}\) contributes \(-\nu(\pi M/L)^4\). An explicit method would require \(\Delta t \sim 1/\vert\mathbf{k}\vert_{\max}^4 \sim L^4/(\nu M^4)\) for stability. On a modest \(64 \times 64\) grid, that's around \(10^{-6}\). We would need a million steps to simulate one time unit. A bit pricey.

The solution: treat the linear part exactly. Write the PDE in Fourier space as

\[\partial_t \hat{u} = \hat{L} \hat{u} + \hat{N}(u),\]

where \(\hat{L} = \text{diag}(\vert\mathbf{k}\vert^2 - \nu\vert\mathbf{k}\vert^4)\). The exact solution over a time step \(h\) satisfies

\[\hat{u}(t+h) = e^{\hat{L}h}\,\hat{u}(t) + \int_0^h e^{\hat{L}(h-\tau)}\,\hat{N}(u(t+\tau))\,d\tau.\]

The first term handles the linear part exactly with no stability restriction from the \(\vert\mathbf{k}\vert^4\) term. The integral involving the nonlinearity is approximated by a Runge–Kutta scheme, giving a family of Exponential Time Differencing methods. The specific one used here is ETDRK4 ⁴, a fourth-order scheme. Its coefficients involve the \(\varphi\)-functions, which arise naturally from the integral in the exact solution formula:

\[\begin{aligned} \varphi_0(z) &= e^z, & \varphi_1(z) &= \frac{e^z - 1}{z}, \\[4pt] \varphi_2(z) &= \frac{e^z - 1 - z}{z^2}, & \varphi_3(z) &= \frac{e^z - 1 - z - z^2/2}{z^3}. \end{aligned}\]

The ETDRK4 step reads:

\[\begin{aligned} a &= \varphi_0(\tfrac{h}{2}\hat{L})\,u_n + \tfrac{h}{2}\varphi_1(\tfrac{h}{2}\hat{L})\,N(u_n), \\ b &= \varphi_0(\tfrac{h}{2}\hat{L})\,u_n + \tfrac{h}{2}\varphi_1(\tfrac{h}{2}\hat{L})\,N(a), \\ c &= \varphi_0(\tfrac{h}{2}\hat{L})\,a + \tfrac{h}{2}\varphi_1(\tfrac{h}{2}\hat{L})\,(2N(b) - N(u_n)), \\ \end{aligned}\] \[u_{n+1} = \varphi_0(h\hat{L}) u_n + h(\varphi_1-3\varphi_2+4\varphi_3)N(u_n) + h(2\varphi_2-4\varphi_3)(N(a)+N(b)) + h(-\varphi_2+4\varphi_3)N(c),\]

where all \(\varphi_{\ell}\) are evaluated at \(h\hat{L}\) in the final line.

Evaluating \(\varphi_1, \varphi_2, \varphi_3\) near \(z = 0\) requires care — the formulas above suffer from catastrophic cancellation. The simulator uses a numerically stable evaluation ⁵, falling back to the explicit formula for \(\vert z\vert > 5\) where it is safe.

The simulator also offers a simpler IMEX BDF2 integrator (implicit-explicit Backward Differentiation Formula, second order). It treats the stiff linear terms implicitly and the nonlinearity explicitly, trading some accuracy for better stability on coarse grids. This is what Reusser's notebook uses.

Why does it still blow up sometimes?

If you push the domain size too large or the hyperviscosity too low, the simulation will eventually blow up and the solver will stop. This is frustrating, especially if you've seen other KS simulators (notably Ricky Reusser's beautiful WebGL notebook) that seem to run forever without a hiccup. What gives?

Half-precision arithmetic as accidental hyperviscosity. Reusser's simulation runs on the GPU using WebGL textures, which default to 16-bit floating point (float16). Every multiplication discards information at high wavenumbers. This acts as an unintentional low-pass filter that continuously damps the very modes most prone to instability. Our solver uses Float64Array (double precision) in a Web Worker, preserving energy across all resolved scales, including the ones that can run away.

A different form of the equation. Reusser follows Kalogirou et al. ⁶, who shift the solution by a constant drift: \(V = u - ct\). This introduces a term proportional to \(c\,dt\) into the BDF2 update formula that provides additional damping of previous timesteps. Our implementation solves the KS equation directly without this shift.

Energy-leaking post-processing. After every timestep, Reusser's code runs an extra cycle to remove imaginary leakage from the complex FFT of real data. Each such pass removes a tiny amount of energy. Our code preserves Hermitian symmetry analytically (if \(\hat{u}_{-\mathbf{k}} = \overline{\hat{u}_{\mathbf{k}}}\), the IFFT is real), so no stripping is needed (but we also don't get the free dissipation).

The IMEX BDF2 integrator mentioned above follows the same approach as Reusser's notebook but without the half-float damping or the \(c\)-shift.

---

Try it yourself

The simulator is live at francesco.ballerin.it/ks_equation. Here are some parameter combinations to get you started:

Regime	\(L\)	\(\nu\)	\(dt\)
Stable	2	1.0	0.25
Cellular	10	1.0	0.25
Chaotic	30	1.0	0.25
Extreme (will rapidly blow up)	30	0.1	0.1

---

References