Goal: derive two derivative approximations and explain why accuracy depends on resolution.
Start with an answer we know
On the periodic interval \([0,2\pi)\), take
\[ u(x)=1+0.1\cos(mx),\qquad u'(x)=-0.1m\sin(mx). \]
Sample at \(x_j=jh\), with \(h=2\pi/N\). Do not repeat the endpoint: \(x_N\) is the same periodic point as \(x_0\). We will compare numerical derivatives at these nodes with the analytic derivative above.
Finite differences: use neighbouring values
Taylor expansion gives
\[ u(x+h)-u(x-h)=2hu'(x)+\tfrac{h^3}{3}u'''(x)+O(h^5). \]
Dividing by \(2h\) gives the centred second-order stencil. Combining neighbours at distances \(h\) and \(2h\) cancels its leading error and gives the fourth-order stencil:
\[ (D_2u)_j=\frac{u_{j+1}-u_{j-1}}{2h},\qquad (D_4u)_j=\frac{u_{j-2}-8u_{j-1}+8u_{j+1}-u_{j+2}}{12h}. \]
Indices wrap periodically. For a fixed sufficiently smooth function, the truncation errors scale like \(h^2\) and \(h^4\). Halving \(h\) should eventually reduce them by factors near 4 and 16, until another error source dominates. Fourth order needs two neighbours on each side; second order needs one. This difference becomes a halo-width requirement when the grid is distributed over MPI ranks.
Fourier differentiation
Represent the periodic samples by a trigonometric interpolant,
\[ u_N(x)=\sum_{k\in K_N}\widehat u_k e^{ikx},\qquad \widehat u_k=\frac1N\sum_{j=0}^{N-1}u_j e^{-ikx_j}. \]
Each basis function has a known derivative:
\[ u_N'(x)=\sum_{k\in K_N}ik\widehat u_k e^{ikx}. \]
Transform, multiply by \(ik\), and transform back. The FFT is an efficient way to calculate the transforms; the differentiation rule comes from the basis. For even \(N\), real first differentiation sets the Nyquist multiplier to zero: the derivative of its cosine interpolant is zero at the sample nodes.
For our single cosine, Fourier differentiation is exact up to roundoff when its mode is resolved below Nyquist. More generally, the rate at which spectral error falls depends on smoothness. Analytic periodic functions can converge geometrically as \(N\) grows. Discontinuities and a mismatch at the periodic boundary remove that advantage; a global Fourier representation can show Gibbs oscillations near a jump.
Experiment with the same samples
Predict what happens when you double \(N\). Then raise \(m\) while keeping \(N\) fixed. Compare errors at the nodes; the lines joining numerical values are only a visual guide.
Interactive experiment / periodic 1D differentiation
Mode m = 3 is below Nyquist. Increase N: FD errors fall while the resolved Fourier derivative approaches roundoff. Interactive controls appear when the experiment is ready.
Browser reference calculation, not OpenPFC solver output. Error is the maximum absolute error at the grid nodes. Drag to zoom; double-click the plot to reset.
Sampling and aliasing
At the nodes, \(e^{i(k+qN)x_j}=e^{ikx_j}\) for any integer \(q\). Samples cannot distinguish these frequencies. A high-frequency signal can therefore masquerade as a lower-frequency one. Increasing the accuracy of a derivative formula cannot recover information absent from its input.
For this experiment require \(m<N/2\). With \(N=16,m=13\), mode 13 aliases to \(-3\) and has the same cosine samples as mode 3. The derivative of that interpolant is \(-0.3\sin(3x)\); the original derivative is \(-1.3\sin(13x)\). At the nodes, \(\sin(13x_j)=-\sin(3x_j)\), so their maximum difference is \(1.6\). This is the sampling checkpoint expressed as a derivative error. At Nyquist, this particular cosine’s true first derivative vanishes at the nodes too; agreement there does not identify the original continuous signal or recover an arbitrary phase.
Nonlinear products can create unresolved high frequencies even when their input fields were resolved. That is why pseudospectral nonlinear solvers need a deliberate aliasing strategy. A fine-looking Fourier curve alone does not establish that a nonlinear calculation is resolved.
Explain the error with a modified wavenumber
Apply each stencil to a single mode \(e^{ikx}\) and write its derivative as \(i\widetilde{k}e^{ikx}\). With \(\theta=kh\),
\[ \widetilde{k}_2=\frac{\sin\theta}{h},\qquad \widetilde{k}_4=\frac{8\sin\theta-\sin(2\theta)}{6h}. \]
For small \(\theta\), these approach \(k\). Near the grid limit they depart from it, and both vanish at \(\theta=\pi\). This explains why many points per wavelength matter even before complete aliasing occurs. Fourier differentiation uses the represented wavenumber directly, subject to the real Nyquist convention above.
Try these checks:
- Hold \(m=1\) and compare \(N=64\) with \(128\). Do the FD error ratios approach 4 and 16? Does the Fourier error remain near roundoff?
- Set \(N=16,m=13\). Explain the Fourier error using the aliased frequency, then increase \(N\) until this mode is resolved.
- Imagine the input is noisy. Which high frequencies does each formula amplify, and why is a small truncation error not the only concern?
Checkpoint: predict the refinement ratios
Hold m = 1 and double N from 64 to 128. Which errors should fall, and by how much?
For fixed smooth input, FD2 error is \(O(h^2)\) and FD4 is \(O(h^4)\). Halving \(h\) gives ratios approaching 4 and 16. The resolved single-mode Fourier error remains near floating-point roundoff rather than following that algebraic trend. Increasing \(m\) along with \(N\) would change the function and would not be this same refinement test.
Connect the idea to OpenPFC
OpenPFC’s finite-difference fields use halos to make neighbouring values available on each rank. Spectral stacks use HeFFTe-backed distributed FFTs. The browser uses a small direct discrete Fourier transform so the reference calculation is transparent; it is not an FFT implementation or a speed comparison.
Use the examples catalogue for runnable C++ examples, the spectral-stack explanation for the transform workflow, and Trefethen’s spectral-methods text for a deeper treatment of smoothness, interpolation and differentiation.