The Spectral Theorem says that every symmetric (or Hermitian) operator can be diagonalized in an orthonormal eigenbasis. For the Laplacian \(-d^2/dx^2\), that eigenbasis is the set of sinusoids. Fourier analysis is therefore just linear algebra in an infinite-dimensional space: expressing a function as a Fourier series is the same operation as expanding a vector in an eigenbasis. This viewpoint makes Fourier analysis less magical and reveals exactly why it is so powerful: it diagonalizes the Laplacian, turning differential equations into simple algebra in frequency space. The Discrete Fourier Transform (DFT) is the finite-dimensional version of this, represented as a unitary matrix. Computing it efficiently via the Fast Fourier Transform (FFT) is one of the most practically impactful algorithms ever discovered.
The Fourier Basis
In Section~above, we found that the eigenfunctions of \(-d^2/dx^2\) on \([0,1]\) are \(u_k(x) = \sqrt{2}\sin(k\pi x)\) with eigenvalues \(k^2\pi^2\). These functions are orthonormal in \(L^2[0,1]\), and they span the whole space. Expanding any function \(f\) in this basis is the Fourier sine series. On the periodic domain \([0, 2\pi]\), the natural eigenbasis is the complex exponentials \(e^{ikx}\), which give the full (complex) Fourier series. We work with the periodic setting here because it connects most directly to the DFT.
For \(f \in L^2([0, 2\pi])\) (periodic), the Fourier series of \(f\) is \[f(x) = \sum_{k=-\infty}^{\infty} \hat{f}_k\, e^{ikx},\] where the Fourier coefficients are \[\hat{f}_k = \frac{1}{2\pi}\int_0^{2\pi} f(x)\,e^{-ikx}\,dx.\] The coefficient \(\hat{f}_k\) is the inner product of \(f\) with the basis function \(e^{ikx}/(2\pi)\), exactly as in a standard change of basis.
The functions \(\phi_k(x) = \frac{1}{\sqrt{2\pi}}e^{ikx}\), \(k \in \fZ\), form an orthonormal set in \(L^2([0, 2\pi])\): \[\langle \phi_j, \phi_k \rangle
= \frac{1}{2\pi}\int_0^{2\pi} e^{ijx}\overline{e^{ikx}}\,dx
= \frac{1}{2\pi}\int_0^{2\pi} e^{i(j-k)x}\,dx
= \begin{cases} 1 & j = k, \\ 0 & j \neq k. \end{cases}\] Moreover, they form a complete orthonormal basis: every \(f \in L^2([0,2\pi])\) can be recovered exactly from its Fourier coefficients via the Fourier series.
Orthonormality. For \(j \neq k\): \[\frac{1}{2\pi}\int_0^{2\pi} e^{i(j-k)x}\,dx
= \frac{1}{2\pi} \cdot \frac{e^{i(j-k)x}}{i(j-k)}\Bigg|_0^{2\pi}
= \frac{e^{2\pi i(j-k)} - 1}{2\pi i(j-k)} = 0,\] since \(e^{2\pi i m} = 1\) for any integer \(m\). For \(j = k\), the integrand is 1 and the integral equals \(2\pi\), giving norm 1 after dividing by \(2\pi\).
Completeness. Showing that every \(L^2\) function is approximated arbitrarily well by finite Fourier sums is a deeper result (the Riesz-Fischer theorem) that goes beyond this course. It is the infinite-dimensional analog of the statement that an orthonormal set in \(\fR^n\) that spans \(\fR^n\) is a basis.
The change-of-basis picture. In \(\fR^n\), changing from the standard basis \(\{\be_j\}\) to an orthonormal basis \(\{\bq_j\}\) is accomplished by the unitary matrix \(\bQ\): the \(j\)-th coordinate of \(\bx\) in the new basis is \(\bq_j^T\bx\). The Fourier transform does exactly the same thing in \(L^2\): the \(k\)-th coordinate of \(f\) in the Fourier basis is \(\hat{f}_k =
\langle \phi_k, f\rangle\). Just as \(\bQ\) diagonalizes a symmetric matrix (\(\bA = \bQ\bD\bQ^T\)), the Fourier transform diagonalizes the Laplacian: \(-\widehat{f''(x)}_k = k^2\hat{f}_k\).
Refer to the result above and the result above.
Compute the Fourier coefficients \(\hat{f}_k\) for \(f(x) = x\) on \([0, 2\pi]\) by evaluating the integral directly. For which \(k\) are the coefficients nonzero?
Verify the diagonalization claim: if \(f(x) = e^{ikx}\), compute \(-f''(x)\) and show that \(-\widehat{f''}_k = k^2\hat{f}_k\) without any integration.
Why does the Fourier transform turn the ODE \(-u'' = f\) into the algebraic equation \(k^2\hat{u}_k = \hat{f}_k\)? When is this equation not solvable?
Circulant Matrices and the Convolution Theorem
The DFT is more than just a change of basis for signals: it simultaneously diagonalizes an entire class of matrices. Circulant matrices arise naturally whenever a system has translation symmetry (same local rule everywhere), which includes the periodic discrete Laplacian and discrete convolution operators. Diagonalizing them with the DFT reduces matrix-vector multiplication from \(O(n^2)\) to \(O(n \log n)\).
A matrix \(\mathbf{C} \in \fC^{n \times n}\) is circulant if each row is a cyclic right-shift of the previous row. It is completely determined by its first column \(\bc = (c_0, c_1, ..., c_{n-1})^T\): \[\mathbf{C} =
\begin{pmatrix}
c_0 & c_{n-1} & c_{n-2} & \cdots & c_1 \\
c_1 & c_0 & c_{n-1} & \cdots & c_2 \\
c_2 & c_1 & c_0 & \cdots & c_3 \\
\vdots & & \ddots & & \vdots \\
c_{n-1} & c_{n-2} & \cdots & c_1 & c_0
\end{pmatrix}.\] The product \(\mathbf{C}\bx\) computes the circular convolution \(\bc * \bx\).
Every circulant matrix \(\mathbf{C}\) with first column \(\bc\) satisfies \[\mathbf{C} = \mathbf{F}^*\,\mathrm{diag}(\hat{\bc})\,\mathbf{F},\] where \(\hat{\bc} = \mathbf{F}\bc\) is the DFT of \(\bc\). The eigenvalues of \(\mathbf{C}\) are exactly the entries of \(\hat{\bc}\) (the DFT coefficients of the first column), and the eigenvectors are the columns of \(\mathbf{F}^*\).
For any \(\bx, \by \in \fC^n\), the circular convolution satisfies \[\mathbf{F}(\bx * \by) = (\mathbf{F}\bx) \odot (\mathbf{F}\by),\] where \(\odot\) denotes entry-wise (Hadamard) product. Circular convolution in the time domain is pointwise multiplication in the frequency domain.
Step 1: Express convolution as a matrix-vector product. Circular convolution \(\bz = \bx * \by\) means \(z_j = \sum_{k=0}^{n-1} x_k y_{(j-k)\bmod n}\). This is exactly \(\bz = \mathbf{C}_{\bx}\by\) where \(\mathbf{C}_{\bx}\) is the circulant matrix with first column \(\bx\).
Step 2: Diagonalize. By the result above, \(\mathbf{C}_{\bx} = \mathbf{F}^*\mathrm{diag}(\hat{\bx})\mathbf{F}\), so: \[\bz = \mathbf{C}_{\bx}\by = \mathbf{F}^*\mathrm{diag}(\hat{\bx})\mathbf{F}\by
= \mathbf{F}^*(\hat{\bx} \odot \hat{\by}).\]
\[\hat{\bz} = \mathbf{F}\bz = \mathbf{F}\mathbf{F}^*(\hat{\bx} \odot \hat{\by})
= \hat{\bx} \odot \hat{\by},\] using \(\mathbf{F}\mathbf{F}^* = \bI\).
Application to the Poisson equation. The 1-D discrete Laplacian with periodic boundary conditions is circulant. Its eigenvalues are the DFT of the first column \(\bc = \frac{1}{h^2}(2, -1, 0, ..., 0, -1)^T\). This gives an \(O(n \log n)\) solver for the periodic Poisson equation:
Compute \(\hat{\mathbf{f}} = \mathbf{F}\mathbf{f}\) via FFT.
Divide entry-wise: \(\hat{u}_k = \hat{f}_k / \hat{c}_k\).
Invert: \(\mathbf{u} = \mathbf{F}^*\hat{\mathbf{u}}\) via inverse FFT.
Compare this to the \(O(n^3)\) dense solve or \(O(n^{3/2})\) sparse solve. On a 2-D grid of size \(n \times n\) (\(N = n^2\) unknowns), the FFT-based solver costs \(O(N \log N)\), making it the method of choice for large periodic problems.
Refer to the result above and the result above.
Build the \(5 \times 5\) circulant matrix with first column \(\bc = (2, -1, 0, 0, -1)^T\) (the periodic 1-D Laplacian for \(n = 5\)). Compute its eigenvalues directly and via \(\hat{\bc} = \mathbf{F}\bc\) using np.fft.fft(c). Verify they match.
Implement circular convolution two ways: (a) directly as a matrix-vector product using scipy.linalg.circulant, and (b) using np.fft.fft and pointwise multiplication. Verify they agree for random \(\bx, \by \in \fR^{16}\).
Implement the FFT-based Poisson solver described in the note above for \(n = 100\) with \(f(x) = \sin(2\pi x)\). The exact solution is \(u(x) = \frac{1}{4\pi^2}\sin(2\pi x)\). Verify your solver recovers it.