Fourier spectral methods and Fourier neural operators

Before this lesson: circulant matrices and the DFT (linear algebra 19); the method of lines, RK4 and its stability region (lesson 02). Complex exponentials: $e^{ikx} = \cos kx + i\sin kx$. Time: about 100 minutes, plus the notebook.

---

Why this matters #

• Weather and climate models. ECMWF's IFS is a spectral transform model: it computes derivatives in spherical-harmonic space (lesson 09 is the spherical version of this lesson). NeuralGCM couples a spectral dynamical core with learned physics.
• Neural operators. The FNO (Li et al., 2021) learns maps between functions, for example from an initial condition to the solution at a later time. FourCastNet's adaptive FNO (AFNO) and the spherical FNO (SFNO) are its descendants in weather models.
• Training data. Many datasets for neural PDE solvers (Burgers, Navier–Stokes, Kuramoto–Sivashinsky) are produced by pseudo-spectral solvers like the one in this lesson.
• One shared trap. A product in a solver and a pointwise activation in a network both create frequencies that the grid cannot hold. Understanding aliasing tells you when "resolution-invariant" is true and when it is not.

1. Predict first #

2. The picture #

Derivatives become multiplications #

A Fourier mode is an eigenfunction of the derivative:

So in Fourier space, $d/dx$ is multiplication by $ik$, and $d^2/dx^2$ is multiplication by $-k^2$. A linear PDE with constant coefficients splits into one scalar ODE per mode. Take the heat equation $u_t = \nu u_{xx}$ with viscosity $\nu > 0$. The coefficient $\hat u_k$ of mode $k$ obeys $\hat u_k' = -\nu k^2 \hat u_k$, so $\hat u_k(t) = e^{-\nu k^2 t}\,\hat u_k(0)$. High modes decay fastest.

On a grid this gives the spectral derivative: FFT, multiply mode $k$ by $ik$, inverse FFT. It uses all $n$ grid points at once. A finite difference uses only neighbours.

Interactive figure: loading…

Spectral accuracy, in words. A smooth periodic function has Fourier coefficients that decay very fast. $n$ grid points hold the lowest $n$ modes, so the error is about the size of the first mode that is left out. For an analytic function the coefficients fall like $e^{-a\lvert k\rvert}$, so the error falls like $e^{-an/2}$: each doubling of $n$ roughly squares the error. For $\sin x$, every grid with $n \ge 3$ gives the exact derivative.

Aliasing: a fast wave can hide as a slow one #

$n$ grid points can tell apart only $n$ modes. On the grid $x_j = 2\pi j/n$, the waves $e^{ikx}$ and $e^{i(k+n)x}$ take the same values at every grid point, because $e^{inx_j} = e^{2\pi i j} = 1$.

Interactive figure: loading…

Products create high modes. A product $u \cdot u$ of two fields with modes $\lvert k\rvert \le K$ has modes up to $2K$. In Fourier space a product is a convolution of coefficient vectors, which costs $O(n^2)$. The pseudo-spectral method instead goes to the grid (inverse FFT), multiplies point by point ($O(n)$) and comes back (FFT): $O(n\log n)$ in total. But the new modes above $n/2$ do not fit on the grid. They fold back onto lower modes and add wrong energy there.

The FNO picture #

A neural operator learns a map between functions, for example from $u(\cdot, 0)$ to $u(\cdot, T)$. The FNO does it in three stages. Lift the input function to $d_v$ channels with a pointwise linear map. Then repeat a Fourier layer: FFT, keep the lowest $k_{\max}$ modes, multiply each kept mode by a learned complex $d_v \times d_v$ matrix, inverse FFT, add a pointwise linear map of the input, and apply a pointwise activation. Finally, project to the output with a pointwise network. The weights belong to wavenumbers, not to grid points. So one set of weights serves every grid. Section 3 says exactly when this gives the same output on every grid.

3. Definitions and main results #

Fourier series on a grid #

Take the domain $[0, 2\pi)$. (For a period $L$, mode $m$ has the wavenumber $k = 2\pi m/L$.) The grid is $x_j = 2\pi j/n$, $j = 0, \dots, n - 1$. A function has the Fourier series $u(x) = \sum_k c_k e^{ikx}$ with coefficients $c_k$. The discrete Fourier transform (DFT), in NumPy's convention, is

For even $n$, the FFT entries come in the order $k = 0, 1, \dots, n/2 - 1, -n/2, \dots, -1$ (np.fft.fftfreq(n, d=1/n) returns exactly these integers). If $u$ has only modes $\lvert k\rvert < n/2$, then $\hat u_k = n\,c_k$.

The spectral derivative of order $p$ is $\operatorname{ifft}\big((ik)^p\,\hat u\big)$. For odd $p$ and even $n$, set the Nyquist mode $k = n/2$ to zero. Reason: the real interpolant holds this mode as $\cos(nx/2)$. Its odd derivatives are multiples of $\sin(nx/2)$, which is zero at every grid point.

Spectral accuracy (a standard result; see Boyd, ch. 2). If $u$ is periodic, analytic and bounded in the strip $\lvert\operatorname{Im} x\rvert < a$, then $\lvert c_k\rvert \le C e^{-a\lvert k\rvert}$, and the error of the spectral derivative on $n$ points falls like $e^{-an/2}$, up to slowly varying factors. If $u$ has only finitely many continuous derivatives, the error falls only like a power of $1/n$. Example: $u = 1/(2 - \cos x)$ has $c_k = (2 - \sqrt3)^{\lvert k\rvert}/\sqrt3 = e^{-a\lvert k\rvert}/\sqrt3$ with $a = \operatorname{arccosh} 2 \approx 1.317$.

Aliasing #

Theorem 1 (aliasing). On the grid $x_j = 2\pi j/n$, $e^{i(k + qn)x_j} = e^{ikx_j}$ for every integer $q$. So if $u(x) = \sum_p c_p e^{ipx}$ (absolutely convergent), the DFT of the samples is

Proof of Theorem 1

1. $e^{iqnx_j} = e^{2\pi i qj} = 1$, because $qj$ is an integer.
2. $u_j = \sum_p c_p e^{ipx_j}$. Group the modes $p$ by their remainder $k = p \bmod n$. By step 1, $e^{ipx_j} = e^{ikx_j}$, so $u_j = \sum_{k=0}^{n-1}\big(\sum_q c_{k+qn}\big)e^{ikx_j}$.
3. The DFT on $n$ points is invertible, and $u_j = \frac1n\sum_k \hat u_k e^{ikx_j}$. Compare coefficients: $\hat u_k/n = \sum_q c_{k+qn}$. $\blacksquare$

So the grid cannot tell mode $k$ from $k \pm n$, $k \pm 2n$, and so on. Energy above the Nyquist mode $n/2$ does not disappear; it shows up as wrong energy at a lower mode.

Theorem 2 (Orszag's 2/3 rule). Let $u$ and $w$ have only modes $\lvert k\rvert \le K$. Multiply them point by point on $n$ grid points, take the FFT, and keep the modes $\lvert k\rvert \le K$. If $3K < n$, the kept coefficients equal those of the exact product $uw$. If $3K \ge n$ (and $K < n$), aliasing can corrupt kept modes.

Why this is true, in words. The product has modes up to $2K$. A mode $p$ above $n/2$ lands at $p - n$. The highest one, $2K$, lands at $2K - n$. If this is still below $-K$, every wrong mode lands outside the kept band, and the truncation throws it away.

Proof of Theorem 2

1. The exact product has modes $\lvert p\rvert \le 2K$; its coefficients are the convolution $\sum_q c^u_q c^w_{p-q}$.
2. Assume $3K < n$, so $2K < n$. A mode with $\lvert p\rvert < n/2$ appears on the grid as itself. A mode with $n/2 \le p \le 2K$ appears at $p - n$ (Theorem 1), and $p - n$ lies in $[-n/2,\ 2K - n]$. Negative $p$ are the mirror image.
3. The largest alias is $2K - n$. It misses every kept mode when $2K - n < -K$, that is, $3K < n$. Then all corrupted entries lie at $\lvert k\rvert > K$, and the truncation removes them.
4. Sharpness: if $3K \ge n$ (and $K < n$), take $u = w = \cos Kx$. The product has mode $p = 2K$, which aliases to $2K - n \in [-K, K)$, a kept mode. $\blacksquare$

In a solver: keep the state band-limited with the mask $\lvert k\rvert < n/3$. Go to the grid, multiply, come back, and apply the mask again. The other classic option is padding (the 3/2 rule): compute the product on a grid with $3/2$ times more points (zero-padded coefficients) and truncate back. It keeps all $n$ modes of the state, at the price of larger FFTs.

Why aliasing matters: energy #

Theorem 3 (energy conservation). Inviscid Burgers $u_t + (u^2/2)_x = 0$ on a periodic domain conserves the energy $E = \frac12\int u^2\,dx$ for smooth solutions. Its Galerkin truncation $u_t = -P_K\,(u^2/2)_x$, where $P_K$ keeps the modes $\lvert k\rvert \le K$ and $u = P_K u$, also conserves $E$.

Proof of Theorem 3

1. $\dfrac{dE}{dt} = \int u\,u_t\,dx = -\int u\,P_K\big(u^2/2\big)_x\,dx$.
2. $P_K$ is an orthogonal projection, so it is self-adjoint: $\int u\,P_K g\,dx = \int (P_K u)\,g\,dx = \int u\,g\,dx$. Hence $\dfrac{dE}{dt} = -\int u\,(u^2/2)_x\,dx = -\int u^2 u_x\,dx$.
3. $u^2 u_x = (u^3/3)_x$, and the integral of the derivative of a periodic function is zero. (Without truncation, step 2 is not needed.) $\blacksquare$

With the 2/3 rule, the pseudo-spectral right-hand side is the Galerkin one (Theorem 2). So the dealiased semi-discrete scheme (before time stepping; RK4 adds a small time-step error) conserves energy for $\nu = 0$. For $\nu > 0$ it gives $dE/dt = -\nu\int u_x^2\,dx \le 0$, as the PDE does. With aliasing, the computed term is $P_K(u^2/2)_x$ plus an aliasing error, step 2 fails, and energy can be created. In an under-resolved run this feeds a pile-up of energy at the highest modes (spectral blocking) and then a blow-up. By Parseval, $E = \frac{\pi}{n^2}\sum_k \lvert\hat u_k\rvert^2$ on the grid, so the energy is easy to monitor.

The pseudo-spectral Burgers solver #

The viscous Burgers equation $u_t + (u^2/2)_x = \nu u_{xx}$ on $[0, 2\pi)$ becomes, mode by mode (method of lines),

with $\widehat{(u^2)}$ computed on the grid and masked as in Theorem 2. RK4 then integrates this system of ODEs (lesson 02).

An exact solution to test against (Cole–Hopf). The substitution $u = -2\nu\,\phi_x/\phi$ turns Burgers into the heat equation $\phi_t = \nu\phi_{xx}$ (derivation in section 8, question 3). For $u(x, 0) = \sin x$, choose $\phi(x, 0) = e^{\cos x/(2\nu)}$; then $-2\nu\phi_x/\phi = \sin x$. With the modified Bessel functions $I_k$ and $R = 1/(2\nu)$, $e^{R\cos x} = I_0(R) + 2\sum_{k\ge1} I_k(R)\cos kx$. Each mode decays as in the heat equation:

$\phi$ stays positive (the heat equation keeps positive data positive), so $u$ has no poles. For small $\nu$ the series loses digits to cancellation; the widget below therefore uses the equivalent heat-kernel integral, which is stable for every $\nu$.

The time step. The viscous term has eigenvalues $-\nu k^2$ (real, negative). The advection term behaves roughly like eigenvalues $\pm i\,k\max\lvert u\rvert$ (imaginary). The RK4 stability region contains $[-2.785, 0]$ on the real axis and $[-2.83i, 2.83i]$ on the imaginary axis (lesson 02). With $K$ the largest kept mode, this gives $h \le 2.785/(\nu K^2)$ and roughly $h \lesssim 2.83/(K\max\lvert u\rvert)$. The viscous limit shrinks like $1/n^2$: the problem is stiff. An integrating factor or ETDRK4 removes this limit (lesson 02).

Interactive figure: loading…

The Fourier neural operator #

Operator learning. We want a map $\mathcal G: a \mapsto u$ between function spaces, for example $u(\cdot, 0) \mapsto u(\cdot, T)$ for Burgers. The data are pairs $(a_i, u_i)$, sampled on grids. Two well-known architectures are DeepONet (Lu et al., 2021), which combines a "branch" network for the input function with a "trunk" network for the output location, and the FNO. Kovachki, Lanthaler and Stuart (2024) review the theory.

The Fourier layer. Let $v_\ell(x) \in \R^{d_v}$ be the $d_v$ channels at the point $x$. Then

with $W \in \R^{d_v \times d_v}$, $b \in \R^{d_v}$, a pointwise activation $\sigma$ (GELU) and a learned matrix-valued kernel $\kappa$. By the convolution theorem, $\widehat{(\mathcal K v)}_k = \hat\kappa_k\,\hat v_k$. The FNO learns $R_k = \hat\kappa_k \in \C^{d_v \times d_v}$ directly for the lowest modes $k < k_{\max}$ and sets all higher modes to zero. On a grid: rfft along $x$, multiply mode $k$ by $R_k$, irfft back to $n$ points. For one channel, $\mathcal K$ is a circulant matrix (linear algebra 19): a global convolution that is also a low-pass filter.

• Parameters: $k_{\max}d_v^2$ complex numbers per layer in 1-D, independent of $n$.
• Cost: $O(d_v\,n\log n)$ for the FFTs plus $O(k_{\max}d_v^2)$ for the mode mixing.
• Shift equivariance: every part commutes with cyclic shifts of the grid. Appending the grid coordinates as extra input channels, as the reference code does, breaks this on purpose.

Theorem 4 (the Fourier layer is discretisation-invariant). Let $v$ be a real function with modes $\lvert k\rvert \le M$ only. On every grid with $n > 2M$ points, $\mathcal K$ maps the samples of $v$ to the samples of the same function $g(x) = \sum_{\lvert k\rvert < k_{\max}} R_k\,c_k\,e^{ikx}$, where $c_k$ are the coefficients of $v$ and $R_{-k} = \overline{R_k}$ (so $g$ is real). We assume $k_{\max} \le n/2$, so the layer never asks for a mode the grid does not have.

Proof of Theorem 4

1. $n > 2M$, so all modes of $v$ lie below the Nyquist mode, and Theorem 1 has no extra terms: $\hat v_k = n\,c_k$ for $\lvert k\rvert \le M$, and $0$ for the other modes.
2. The layer multiplies $\hat v_k$ by $R_k$ for $\lvert k\rvert < k_{\max}$ and zeroes the rest. The $R_k$ do not depend on $n$. Modes that a coarse grid cannot hold have $c_k = 0$ anyway.
3. The inverse DFT divides by $n$: $\text{output}_j = \frac1n\sum_k R_k\,n\,c_k\,e^{ikx_j} = g(x_j)$. The factor $n$ of step 1 cancels. $\blacksquare$

The pointwise parts ($W$, $b$, lifting, projection) act on the values at each point, so they do not depend on the grid either.

Why the full FNO is only approximately invariant. $\sigma$ is applied point by point, and $\sigma(v)$ is not band-limited. By Theorem 1 its samples alias: on a coarse grid the next layer sees $c_k$ plus an aliasing error, on a fine grid almost only $c_k$. So the outputs at two resolutions differ by an error that falls as the grid gets finer, as fast as the Fourier coefficients of $\sigma(v)$ decay. In addition, if the data have content above the coarse Nyquist mode, the coarse samples are aliased before the network sees them. And "zero-shot super-resolution" cannot create real detail above the modes the network learned.

Interactive figure: loading…

Alias-aware neural operators

The remedy for aliasing in a network is the same as in a solver: apply the nonlinearity on a finer grid and filter back. The convolutional neural operator (CNO; Raonić et al., 2023) upsamples before each activation and low-pass filters after it, and Bartolucci et al. (NeurIPS 2023) define "representation-equivalent" operators for which the discrete and the continuous network agree. In image generation, StyleGAN3 uses the same idea to remove texture sticking. The price is larger FFTs, as with padding in a solver.

Answers to section 1. (1) $u^2 = \frac12 + \frac12\cos 6x$. On 8 points, mode 6 appears as $6 - 8 = -2$, so the FFT shows $\frac12 + \frac12\cos 2x$: a false mode 2. (2) The weights do not change: they live on the modes $k < k_{\max}$. Each Fourier layer computes the same function on both grids for band-limited input (Theorem 4). The full network is only approximately the same, because the activation aliases differently on the two grids.

4. Worked example #

(a) Aliasing by hand. $n = 8$, $u = \cos 3x$. The exact square is $\frac12 + \frac12\cos 6x$. At the grid points, $\cos 6x_j = \cos(2\pi \cdot 6j/8) = \cos(2\pi\cdot(-2)j/8) = \cos 2x_j$. The FFT reports $\frac12 + \frac12\cos 2x$, that is, $\hat u_k/n = \frac14$ at $k = \pm 2$. A solver that keeps all modes below the Nyquist mode ($\lvert k\rvert \le 3$) adds this false $\cos 2x$ to the solution. The 2/3 rule keeps only $\lvert k\rvert < 8/3$, that is, $\lvert k\rvert \le 2$: then $\cos 3x$ is never part of the state.

(b) A dealiased product. $u = \cos x + \cos 2x$ (modes $\lvert k\rvert \le 2$, and $3 \cdot 2 = 6 < 8$). Exactly, $u^2 = \cos^2 x + 2\cos x\cos 2x + \cos^2 2x$. With $\cos^2 a = \frac12(1 + \cos 2a)$ and $2\cos a\cos b = \cos(a - b) + \cos(a + b)$:

On 8 points, mode 3 stays mode 3 and mode 4 is the Nyquist mode. Nothing lands on $\lvert k\rvert \le 2$. After the mask: $1 + \cos x + \frac12\cos 2x$, exactly the kept part of the true product.

(c) A spectral derivative by hand. $u = \sin x + \frac12\cos 2x$ on $n = 8$ points. Since $\sin x = (e^{ix} - e^{-ix})/(2i)$, the FFT gives $\hat u_1 = 8/(2i) = -4i$ and $\hat u_{-1} = 4i$. The cosine gives $\hat u_{\pm 2} = 8 \cdot \frac12 \cdot \frac12 = 2$. Multiply by $ik$:

• mode $1$: $i \cdot (-4i) = 4$; mode $-1$: $(-i)\cdot 4i = 4$;
• mode $2$: $2i \cdot 2 = 4i$; mode $-2$: $-2i \cdot 2 = -4i$.

The inverse FFT divides by 8: $(4e^{ix} + 4e^{-ix})/8 = \cos x$ and $(4ie^{2ix} - 4ie^{-2ix})/8 = -\sin 2x$. So $u' = \cos x - \sin 2x$, exact at every grid point.

(d) The time step. $n = 64$, $\nu = 0.1$, $\max\lvert u\rvert \le 1$. The 2/3 rule keeps $\lvert k\rvert \le 21$ ($3 \cdot 21 = 63 < 64$). Viscous limit: $h \le 2.785/(0.1 \cdot 21^2) \approx 0.063$. Advective limit: $h \lesssim 2.83/21 \approx 0.13$. The viscous term decides. Without dealiasing ($K = 31$) the viscous limit would be $2.785/(0.1 \cdot 31^2) \approx 0.029$. Doubling $n$ divides the viscous limit by about 4.

(e) FNO weights. With $k_{\max} = 16$ and $d_v = 32$, one Fourier layer has $16 \cdot 32 \cdot 32 = 16{,}384$ complex weights (plus $32 \cdot 32 + 32 = 1{,}056$ real numbers for $W$ and $b$), on a grid of 64 points or of 4096 points.

5. In code #

NumPy and SciPy check the claims of this lesson. First, aliasing on 8 points (worked examples a and b):

import numpy as np

n = 8
x = 2 * np.pi * np.arange(n) / n                   # periodic grid on [0, 2π)
k = np.fft.fftfreq(n, d=1 / n).astype(int)         # the integer wavenumber of each FFT entry
print(k)                                           # [ 0  1  2  3 -4 -3 -2 -1]
print(np.allclose(np.cos(6 * x), np.cos(2 * x)))   # True: mode 6 looks like mode 6 - 8 = -2

coef = lambda w: (np.fft.fft(w) / n).real.round(3) + 0.0   # real signals that are even in x
print(coef(np.cos(3 * x) ** 2))   # [0.5  0.   0.25 0.   0.   0.   0.25 0.  ]
                                  # exact: 1/2 + cos(6x)/2, so the 0.25 at k = ±2 is aliasing
u = np.cos(x) + np.cos(2 * x)     # modes |k| <= 2: the 2/3-rule band for n = 8
print(coef(u ** 2))               # [1.   0.5  0.25 0.5  0.5  0.5  0.25 0.5 ]
                                  # modes |k| <= 2 (1, 0.5, 0.25) are exact

Second, spectral accuracy against second-order finite differences, for $u = 1/(2 - \cos x)$:

def spectral_derivative(u):
    n = len(u)
    k = np.fft.fftfreq(n, d=1 / n)
    k[n // 2] = 0                                   # Nyquist mode: zero for an odd derivative
    return np.fft.ifft(1j * k * np.fft.fft(u)).real

a = np.arccosh(2.0)                                 # analytic in the strip |Im x| < a
for n in [8, 16, 32, 64]:
    x = 2 * np.pi * np.arange(n) / n
    u, du = 1 / (2 - np.cos(x)), -np.sin(x) / (2 - np.cos(x)) ** 2
    fd = (np.roll(u, -1) - np.roll(u, 1)) / (2 * 2 * np.pi / n)   # periodic central differences
    print(n, f"{np.abs(spectral_derivative(u) - du).max():.1e}",
          f"{np.abs(fd - du).max():.1e}", f"{np.exp(-a * n / 2):.1e}")
# 8 1.3e-02 1.0e-01 5.2e-03
# 16 1.4e-04 4.2e-02 2.7e-05
# 32 7.5e-09 1.1e-02 7.1e-10
# 64 5.0e-15 2.8e-03 5.0e-19

The spectral error follows $e^{-an/2}$ (last column) up to a slowly varying factor, so each doubling of $n$ roughly squares it, until it reaches round-off at $n = 64$. The finite-difference error falls only by about 4 per doubling (order 2).

Third, the pseudo-spectral Burgers solver, checked against Cole–Hopf, and the energy with and without the 2/3 rule:

from scipy.special import ive

def burgers_rhs(uh, k, nu, mask):
    """d(u_hat)/dt for u_t + (u^2/2)_x = nu u_xx, pseudo-spectral, with a mode mask."""
    u = np.fft.ifft(uh).real
    return mask * (-0.5j * k * np.fft.fft(u * u) - nu * k**2 * uh)

def solve_burgers(u0, nu, T, h, dealias=True):
    n = len(u0)
    k = np.fft.fftfreq(n, d=1 / n)                     # integer wavenumbers
    mask = (np.abs(k) < n / 3) if dealias else (np.abs(k) < n / 2)
    uh = mask * np.fft.fft(u0)
    for _ in range(int(round(T / h))):                 # RK4 (lesson 02)
        k1 = burgers_rhs(uh, k, nu, mask)
        k2 = burgers_rhs(uh + h / 2 * k1, k, nu, mask)
        k3 = burgers_rhs(uh + h / 2 * k2, k, nu, mask)
        k4 = burgers_rhs(uh + h * k3, k, nu, mask)
        uh = uh + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    return np.fft.ifft(uh).real

def cole_hopf(x, t, nu, K=200):
    """Exact solution for u(x, 0) = sin x: phi = I0(R) + 2 sum I_k(R) e^{-nu k^2 t} cos kx, R = 1/(2 nu)."""
    k = np.arange(K)
    c = ive(k, 1 / (2 * nu)) * np.exp(-nu * k**2 * t)  # scaled Bessel functions: the scale cancels
    c[1:] *= 2
    phi = c @ np.cos(np.outer(k, x))
    phi_x = -(c * k) @ np.sin(np.outer(k, x))
    return -2 * nu * phi_x / phi

for n in [32, 64, 128]:
    x = 2 * np.pi * np.arange(n) / n
    u = solve_burgers(np.sin(x), nu=0.1, T=1.0, h=1e-3)
    print(n, f"{np.abs(u - cole_hopf(x, 1.0, 0.1)).max():.1e}")
# 32 4.6e-03
# 64 3.7e-05
# 128 2.5e-09

energy = lambda u: 0.5 * np.mean(u**2) * 2 * np.pi   # E = (1/2) ∫ u^2 dx
x = 2 * np.pi * np.arange(64) / 64
for T in [1.0, 1.5]:
    on = solve_burgers(np.sin(x), nu=0.005, T=T, h=2e-3, dealias=True)
    off = solve_burgers(np.sin(x), nu=0.005, T=T, h=2e-3, dealias=False)
    print(T, f"{energy(on):.3f} {energy(off):.3f}")
# 1.0 1.546 1.548
# 1.5 1.409 2.591      (exact: 1.546 at t = 1 and 1.337 at t = 1.5; E(0) = π/2 = 1.571)

Fourth, the two facts behind discretisation invariance: a band-limited input has the same coefficients on every fine-enough grid (Theorem 4), and a pointwise GELU makes coarse grids alias:

from scipy.special import erf

gelu = lambda a: 0.5 * a * (1 + erf(a / np.sqrt(2)))
grid = lambda n: 2 * np.pi * np.arange(n) / n
coeffs = lambda vals, kmax=4: np.fft.rfft(vals)[: kmax + 1] / len(vals)   # c_0 .. c_kmax

v = lambda x: 2 * np.sin(x) + np.cos(3 * x)            # band-limited: modes 1 and 3
print(np.round(coeffs(v(grid(16))), 3) + 0)            # [0. +0.j 0. -1.j 0. +0.j 0.5+0.j 0. +0.j]
print(np.abs(coeffs(v(grid(16))) - coeffs(v(grid(64)))).max() < 1e-14)   # True

ref = coeffs(gelu(v(grid(4096))))                      # practically exact
for n in [16, 32, 64]:
    print(n, f"{np.abs(coeffs(gelu(v(grid(n)))) - ref).max():.1e}")
# 16 4.7e-03
# 32 2.0e-06
# 64 9.0e-15

6. Where it shows up #

• Weather and climate. Orszag published the 2/3 rule in the Journal of the Atmospheric Sciences (1971). Spectral transform models such as ECMWF's IFS, and NeuralGCM's dynamical core, use the spherical version of this lesson (lesson 09). FNO-type emulators (FourCastNet's AFNO, SFNO) learn global dynamics in Fourier or spherical-harmonic space (lesson 16).
• Seismology. Fourier pseudo-spectral methods are a classic tool for seismic wave propagation. Neural operators have been trained to emulate 2-D acoustic wave propagation for any velocity model and source location, and reverse-mode automatic differentiation through them gives the gradients for tomography (Yang et al., The Seismic Record, 2021): the adjoint idea of lesson 04.
• PINN baselines. The Burgers solver of this lesson is the reference against which lesson 11 measures physics-informed networks.
• Deep learning. An FNO layer is a global convolution computed with the FFT; long-convolution sequence models use the same trick. Strided convolutions and pooling in CNNs alias too. Alias-free generators (StyleGAN3) upsample before the nonlinearity and low-pass after it, the same remedy as dealiasing.

7. Common mistakes #

• Wavenumbers in the wrong order, or without the factor $2\pi/L$ for a period $L \ne 2\pi$. Use np.fft.fftfreq(n, d=L/n) * 2 * np.pi.
• A nonzero Nyquist mode in an odd derivative. The result is no longer the derivative of a real interpolant.
• Keeping $\lvert k\rvert \le n/3$ instead of $\lvert k\rvert < n/3$. For $n = 96$ this keeps mode 32, and mode 64 aliases onto $-32$.
• Masking only after the product. The state itself must be band-limited, so mask the initial condition and every right-hand side.
• Blaming the nonlinear term for a blow-up that is really the viscous stability limit $h \le 2.785/(\nu K^2)$: check the time step first, then the energy.
• FNO: inconsistent FFT normalisation (norm="ortho" in one place, the default in another). The output then scales with $n$, and discretisation invariance is lost.
• FNO: asking for more modes than the grid has (rfft gives $n/2 + 1$ modes for even $n$), or evaluating on non-periodic data without padding.

8. Check yourself #

1. On $n = 12$ grid points, which mode does $\cos 10x$ look like?
   Answer

   $10 - 12 = -2$: it looks like $\cos 2x$ (cosine is even, so modes $\pm 2$ give the same real wave).
2. For $n = 48$, what is the largest kept $\lvert k\rvert$ under the 2/3 rule? Show what goes wrong if you keep one more mode.
   Answer

   Keep $\lvert k\rvert < 16$, so $K = 15$ ($3 \cdot 15 = 45 < 48$). With $K = 16$, the product has mode $32 > 24$, which aliases to $32 - 48 = -16$: a kept mode.
3. Derive the Cole–Hopf transform: show that $u = -2\nu\phi_x/\phi$ solves $u_t + uu_x = \nu u_{xx}$ when $\phi_t = \nu\phi_{xx}$.
   Answer

   Write $u = \psi_x$ with $\psi = -2\nu\log\phi$. Burgers is $\partial_x\big(\psi_t + \tfrac12\psi_x^2 - \nu\psi_{xx}\big) = 0$, so it is enough that $\psi_t + \tfrac12\psi_x^2 - \nu\psi_{xx} = 0$. With $\psi_t = -2\nu\phi_t/\phi$, $\psi_x = -2\nu\phi_x/\phi$ and $\psi_{xx} = -2\nu(\phi_{xx}/\phi - \phi_x^2/\phi^2)$: the $\phi_x^2/\phi^2$ terms cancel ($\tfrac12 \cdot 4\nu^2 = 2\nu^2$ against $-\nu \cdot 2\nu$), and the bracket equals $-\tfrac{2\nu}{\phi}(\phi_t - \nu\phi_{xx})$, which is zero for a heat-equation solution. $\phi$ must stay positive, or $u$ has poles.
4. Why is a Fourier layer discretisation-invariant, while a 3-tap convolution kernel learned at 64 points is not?
   Answer

   The Fourier weights $R_k$ are attached to physical wavenumbers $k$, and the DFT scaling cancels (Theorem 4). The 3 weights of a local kernel are attached to the grid spacing: at 256 points they cover a 4 times shorter distance, so they are a different operator with a different frequency response.
5. Why does the dealiased pseudo-spectral scheme for $\nu = 0$ conserve $\sum_k\lvert\hat u_k\rvert^2$, and the aliased one not?
   Answer

   Dealiased, the right-hand side equals the Galerkin term $-P_K(u^2/2)_x$ exactly (Theorem 2). Then $dE/dt = -\int u\,(u^2/2)_x\,dx = -\int (u^3/3)_x\,dx = 0$ (Theorem 3), and $E$ is proportional to $\sum\lvert\hat u_k\rvert^2$ (Parseval). Aliased, the computed term contains an aliasing error whose product with $u$ does not integrate to zero.
6. You double $n$ in the explicit RK4 Burgers solver at fixed $\nu$. What happens to the stable time step, and how can you avoid that?
   Answer

   $K$ doubles, so the viscous limit $2.785/(\nu K^2)$ falls by a factor 4 (the advective limit only by 2). Treat the stiff linear term exactly or implicitly: an integrating factor, ETDRK4 or an IMEX scheme (lesson 02). Then the step is limited mainly by advection.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. How do you compute the first derivative of periodic samples with the FFT?

   Show answer

2. What does spectral accuracy mean?

   Show answer

3. State the aliasing rule on a grid of $n$ points.

   Show answer

4. On $n = 8$ points, where does mode 6 of $u^2$ appear?

   Show answer

5. State Orszag's 2/3 rule.

   Show answer

6. Why does aliasing matter for energy in Burgers?

   Show answer

7. What is the Cole–Hopf transform?

   Show answer

8. Write the FNO Fourier layer.

   Show answer

9. Why is the linear Fourier layer discretisation-invariant?

   Show answer

10. Why is the full FNO only approximately resolution-invariant?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What is aliasing, and how do you avoid it in a pseudo-spectral code?" — Theorem 1, then the 2/3 rule or padding, with the energy argument as the reason it matters.
• "Is an FNO really resolution-invariant?" — Each Fourier layer is, for band-limited inputs (Theorem 4). The pointwise activations alias on coarse grids, so the network is only approximately invariant; test the error at several resolutions.
• "FNO, U-Net or graph network for a PDE surrogate?" — FNO: global receptive field in each layer, $O(n\log n)$, periodic uniform grids, smooth fields (sharp fronts need many modes). U-Net: local and multiscale, better at sharp features, tied to one resolution. Graph or mesh networks: irregular geometry (GraphCast uses an icosahedral mesh). The sphere: SFNO (lesson 09).
• "Why are spectral methods so accurate, and when do they fail?" — Exponential convergence for smooth periodic functions; trouble with shocks (Gibbs), non-periodic boundaries (use Chebyshev points) and complex geometry.
• "Your spectral solver blows up after 100 steps. How do you debug it?" — Check the time step against the viscous limit, check dealiasing, monitor the energy and the spectrum at the highest modes.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It has four exercises, each with a check cell and a folded solution:

1. FFT wavenumbers and the spectral derivative; compare its error with finite differences.
2. Aliasing of a product and the 2/3-rule mask; find the false modes of $\cos^2 3x$ on 8 points.
3. A pseudo-spectral Burgers solver with RK4, checked against Cole–Hopf; energy with and without the mask.
4. Learn a linear Fourier multiplier (the linear part of an FNO layer) from heat-equation data on 64 points by least squares, and test it on 256 points.

11. Watch and read #

The Fourier Transform and DerivativesSteve Brunton · 11 min

Why a derivative becomes multiplication by ik in Fourier space.

Computing Derivatives with FFT [Python]Steve Brunton · 11 min

The spectral derivative in code, compared with finite differences, as in section 2.

Solving PDEs with the FFT [Python]Steve Brunton · 15 min

The method of lines in Fourier space, the step from the derivative to a PDE solver.

Shannon Nyquist Sampling TheoremSteve Brunton · 17 min

Sampling and aliasing from the signal-processing side; compare with Theorem 1.

Fourier Neural Operator (FNO) [Physics Informed Machine Learning]Steve Brunton · 18 min

The FNO architecture and how it is trained, with the Fourier layer of section 3.

ETH Zürich AISE: Fourier Neural OperatorsCAMLab, ETH Zürich · 84 min

A university lecture on FNOs, including aliasing and what resolution invariance really means.

• Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed. (free PDF from the author): §2.2–2.4 and §2.8 (orders of convergence; singularities set the rate); ch. 11 "Aliasing, spectral blocking, and blow-up": §11.2 (aliasing), §11.3 (spectral blocking), §11.5 (the two-thirds rule).
• Driscoll and Braun, Fundamentals of Numerical Computation (free online): §9.5 (trigonometric interpolation), §11.2–11.4 (method of lines, absolute stability, stiffness).
• Li et al., Fourier Neural Operator for Parametric Partial Differential Equations (ICLR 2021): §4 (the Fourier layer), §5.1 (Burgers).
• Kovachki, Lanthaler and Stuart, Operator Learning: Algorithms and Analysis (2024): §1.1 (vectors versus functions), §3.3 (FNO), §4.2 (universal approximation).
• Lu et al., DeepONet (Nature Machine Intelligence, 2021): the branch–trunk operator network.
• Raonić et al., Convolutional Neural Operators (NeurIPS 2023): alias-aware operator layers.
• neuraloperator, theory guide "Fourier Neural Operators": the sections on the Fourier layer and resolution invariance.