Differentiation that costs no accuracy
Back in the finite-difference rung you approximated f'(x) with a central difference like (f(x+h) - f(x-h)) / (2*h). It was honest work, but its error was only O(h^2): halve the spacing and the error drops by a mere factor of four. You could chase higher order with wider stencils, yet each one still buys you a fixed power of h. Spectral methods make a bolder bet. If the function is smooth — infinitely differentiable and, ideally, periodic — its error can fall faster than h^p for every p at once. We call that exponential (or spectral) accuracy, and the engine that delivers it is the FFT.
The idea is one clean trick you already half-know. The DFT writes a sampled periodic function as a sum of waves, each with a known frequency k. Differentiating a wave is trivial: the derivative of e^{i*k*x} is just i*k times e^{i*k*x}. So to differentiate f, you do not touch f in physical space at all. You take its FFT to land in frequency space, multiply every coefficient by i*k, and transform back. Three steps, total cost O(N log N), and the answer is the derivative of the unique trigonometric interpolant through your samples.
spectral derivative of f sampled at N points: F = fft(f) # to frequency space, O(N log N) ik = i * k # k = [0, 1, ..., N/2, -N/2+1, ..., -1] Fp = ik * F # differentiate each wave fp = real(ifft(Fp)) # back to physical space # fp ~ f' to (near) machine precision when f is smooth & periodic
Why smoothness buys exponential accuracy
Here is the heart of it. The error of any Fourier-based method is governed by how fast the function's Fourier coefficients shrink as the frequency k grows — the tail you throw away when you keep only N of them. And that decay rate is dictated by smoothness. If f has a kink (a jump in some derivative), its coefficients fall only like a power of 1/k, so truncation leaves an algebraic O(1/N^p) error — no better than a good finite-difference scheme. But if f is infinitely smooth and periodic, the coefficients decay faster than any power of 1/k, often like e^{-c*N}. Truncating that tail leaves an error that also shrinks exponentially. Smoothness is not a nice-to-have here; it is the whole reason the method is fast.
A second, equally honest reason: a spectral derivative uses every sample point to estimate the derivative at each location, with weights that fall off slowly. A central difference uses only the two nearest neighbours and throws the rest away. By listening to the whole grid, the spectral method extracts far more of the information a smooth function actually carries — which is exactly why it converges so much faster. The flip side is that there is no locality: change one sample and every derivative value shifts a little.
From differentiating to solving
Once you can differentiate by multiplying in frequency space, whole differential equations become almost algebraic. Take the periodic Poisson equation u''(x) = f(x). In physical space it is a differential equation; in Fourier space it turns into i*k applied twice, that is, multiply the coefficient of u by -(k^2). So each Fourier mode satisfies -(k^2) * u_hat_k = f_hat_k, and you simply divide: u_hat_k = -f_hat_k / k^2 (the k=0 mode is fixed by a constant). FFT the right-hand side, divide mode by mode, inverse-FFT, done. A differential operator became diagonal — element-by-element multiplication — which is the cleanest situation in all of numerical computing.
- Sample f on N equally spaced points over one period — the natural grid for a periodic problem.
- Take the FFT to get coefficients f_hat_k for each frequency k.
- Solve mode by mode: u_hat_k = -f_hat_k / k^2 for k not 0, and fix u_hat_0 by the chosen constant.
- Inverse-FFT the u_hat_k back to physical space to read off u — the whole solve is two FFTs and one divide.
This is what a spectral method is: represent the unknown in a basis of global functions (here, Fourier waves), turn the PDE into equations on the coefficients, solve those, and transform back. For time-dependent equations you combine this spatial machinery with an ODE solver from the earlier rung — but beware, the i*k and -(k^2) factors make the system extremely stiff at high resolution, so the very accuracy that makes space cheap can force you toward implicit or specially split time-steppers. Diagonal in space does not mean free in time.
When the problem is not periodic: Chebyshev
Fourier spectral methods assume periodicity, but most real problems live on an interval with boundary conditions and no wrap-around. The naive fix — fit a single high-degree polynomial through equally spaced points — is a trap you have met before: equispaced polynomial interpolation suffers Runge's phenomenon, where the interpolant oscillates wildly near the ends and the error actually grows with degree. So equispaced sampling, which was perfect for periodic Fourier, is poison for non-periodic polynomials.
The cure is to sample at Chebyshev nodes — points clustered near the two ends, the projection onto the x-axis of equally spaced points on a semicircle. On those nodes a high-degree polynomial interpolant is beautifully behaved and, for a smooth function, converges exponentially, just like Fourier did. This is the Chebyshev spectral method, the workhorse for non-periodic domains. And the bridge to today's whole rung is gorgeous: the substitution x = cos(theta) turns a Chebyshev series into a Fourier cosine series, so you can compute Chebyshev coefficients and derivatives with the FFT — the same O(N log N) machine, reused.
The honest fine print
Spectral methods are spectacular but not universal, and a clear-eyed user keeps the trade-offs in view. They demand smooth solutions and simple geometries — a periodic box or an interval, not a wing or an engine block; on complicated shapes the finite-element world from the FEM rung wins on flexibility. They are global, so a local feature touches every coefficient and you cannot easily refine just one region. And the dense, all-to-all coupling can make the linear algebra worse-conditioned than a sparse finite-difference matrix, so watch the condition number just as carefully here as anywhere.
Keep the usual numerical honesty too. Exponential convergence describes how error falls with N until it hits the round-off floor — past that point, adding modes only adds noise, exactly like driving a finite-difference h too small. The reported accuracy is still accuracy = conditioning times stability: a well-conditioned spectral solve loses few digits, an ill-conditioned one loses many regardless of how elegant the transform is. And every number out is an approximation of a real-arithmetic ideal you never actually computed.
Step back and see what the FFT rung has built. You started by writing a signal as a sum of waves with the DFT and its O(N^2) cost; Cooley-Tukey cut that to O(N log N) by divide-and-conquer; aliasing and the Nyquist limit told you which frequencies your samples can honestly represent; the convolution theorem turned a slow O(N^2) sum into a fast O(N log N) one and even multiplied polynomials and big integers for free. This last guide closes the loop: the very same transform that filters audio and compresses JPEGs also differentiates smooth functions to near machine precision and solves their differential equations. One algorithm, an astonishing range — and now it is yours.