Chebyshev Pseudospectral Methods¶

Chebyshev pseudospectral methods extend the spectral approach to non-periodic problems on bounded intervals, achieving the same exponential accuracy as Fourier methods without requiring periodicity.

1. Why Chebyshev?¶

Fourier methods require periodic boundary conditions. For problems on bounded domains with Dirichlet or Neumann boundary conditions, a naive periodic extension introduces discontinuities and destroys spectral accuracy.

Chebyshev methods solve this by using a basis of Chebyshev polynomials that are naturally suited to \([-1, 1]\). They inherit exponential convergence for smooth functions and offer a DCT-based transform with \(O(N \log N)\) complexity.

Key advantages over Fourier for non-periodic problems:

No periodic extension required
Naturally incorporates boundary conditions
Gauss-Lobatto nodes include the endpoints \(\pm 1\)
Spectral accuracy maintained throughout

2. Chebyshev Polynomials¶

The Chebyshev polynomials of the first kind \(T_n : [-1,1] \to [-1,1]\) are defined by the recurrence:

\[T_0(x) = 1, \quad T_1(x) = x, \quad T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)\]

or equivalently, using the trigonometric representation:

\[T_n(x) = \cos(n \arccos x), \quad x \in [-1, 1]\]

The first few polynomials are:

\[T_0 = 1, \quad T_1 = x, \quad T_2 = 2x^2 - 1, \quad T_3 = 4x^3 - 3x, \quad T_4 = 8x^4 - 8x^2 + 1\]

Orthogonality¶

Chebyshev polynomials are orthogonal with respect to the weight \(w(x) = (1-x^2)^{-1/2}\):

\[\int_{-1}^{1} T_m(x) T_n(x) \frac{dx}{\sqrt{1-x^2}} = \begin{cases} 0 & m \neq n \\ \pi/2 & m = n \neq 0 \\ \pi & m = n = 0 \end{cases}\]

3. Chebyshev-Gauss-Lobatto Nodes¶

The Chebyshev-Gauss-Lobatto (CGL) nodes are the extrema of \(T_N\) plus the endpoints:

\[x_j = \cos\!\left(\frac{j\pi}{N}\right), \quad j = 0, 1, \ldots, N\]

These \(N+1\) points lie in \([-1, 1]\) and include both endpoints \(x_0 = +1\) and \(x_N = -1\), which is essential for imposing boundary conditions.

Node Clustering at Boundaries¶

The CGL nodes cluster quadratically near the boundaries:

\[x_j \approx 1 - \frac{j^2\pi^2}{2N^2} \quad (j \text{ small}), \quad \Delta x_{\min} = O(N^{-2})\]

This clustering is the Chebyshev method's mechanism for controlling the Runge phenomenon — the oscillations near boundaries that plague high-degree polynomial interpolation on uniform grids. The boundary clustering means that the condition number of the differentiation matrix grows as \(O(N^2)\) rather than exponentially.

4. The Chebyshev Transform¶

A function \(u\) expanded in Chebyshev polynomials:

\[u(x) = \sum_{n=0}^{N} \hat{u}_n T_n(x)\]

The Chebyshev coefficients can be computed via the Discrete Cosine Transform (DCT). On the CGL nodes:

\[\hat{u}_n = \frac{2}{N c_n} \sum_{j=0}^{N} \frac{u(x_j)}{c_j} \cos\!\left(\frac{n j \pi}{N}\right)\]

where \(c_0 = c_N = 2\), \(c_j = 1\) otherwise. This is precisely the DCT-I and can be computed in \(O(N \log N)\) using FFT.

5. Spectral Differentiation via the Differentiation Matrix¶

The Chebyshev differentiation matrix \(\mathbf{D}\) maps function values at CGL nodes to derivative values:

\[\mathbf{u}' = \mathbf{D} \mathbf{u}\]

where \(\mathbf{u} = [u(x_0), \ldots, u(x_N)]^T\).

The entries of \(\mathbf{D}\) are given by (Trefethen, 2000):

\[D_{jk} = \begin{cases} \displaystyle\frac{c_j}{c_k} \frac{(-1)^{j+k}}{x_j - x_k} & j \neq k \\[6pt] \displaystyle -\frac{x_j}{2(1 - x_j^2)} & 0 < j < N \\[6pt] \displaystyle\frac{2N^2 + 1}{6} & j = k = 0 \\[6pt] \displaystyle-\frac{2N^2 + 1}{6} & j = k = N \end{cases}\]

Higher-order derivatives are computed by matrix powers: \(\mathbf{u}^{(n)} = \mathbf{D}^n \mathbf{u}\).

The condition number of \(\mathbf{D}\) scales as \(O(N^2)\) and of \(\mathbf{D}^n\) as \(O(N^{2n})\).

6. Chebyshev Derivative via Recurrence Relation¶

An alternative to the differentiation matrix uses the recurrence relation on Chebyshev coefficients:

\[c_{n-1} \hat{u}'_{n-1} = \hat{u}'_{n+1} + 2n \hat{u}_n, \quad n = N, N-1, \ldots, 1\]

starting from \(\hat{u}'_N = 0\) and \(\hat{u}'_{N+1} = 0\). This avoids matrix-vector products and runs in \(O(N)\) operations after the initial DCT.

Unrolling the recurrence gives the closed form

\[\hat{u}'_k = \frac{2}{c_k} \sum_{\substack{j > k \\ j + k \text{ odd}}} j\, \hat{u}_j,\]

i.e. a reverse cumulative sum over each parity class. spectraldiffx evaluates it this way (chebyshev_derivative_coeffs), so the step is parallel rather than a sequential scan. Pass method="fft" to ChebyshevDerivative1D/2D/3D to use this path; the default "matrix" path is usually faster for moderate \(N\) because a dense mat-vec maps well onto BLAS.

Full algorithm:

\(\hat{u}_n \leftarrow \text{DCT}(u_j)\)
Apply recurrence to get \(\hat{u}'_n\)
\(u'_j \leftarrow \text{IDCT}(\hat{u}'_n)\)

7. Boundary Conditions¶

The CGL grid includes the endpoints, so boundary conditions are imposed by replacing the first or last row of the differentiation matrix system.

Dirichlet Boundary Conditions¶

For \(u(-1) = a\) and \(u(1) = b\):

Replace the rows for \(j=0\) and \(j=N\) in the system \(\mathbf{D}^2 \mathbf{u} = \mathbf{f}\) with:

\[u_0 = b, \quad u_N = a\]

In matrix form, this sets the first and last rows of the coefficient matrix to \([1, 0, \ldots, 0]\) and \([0, \ldots, 0, 1]\) respectively, with the right-hand side set to \(b\) and \(a\).

Neumann Boundary Conditions¶

For \(u'(-1) = a\) and \(u'(1) = b\):

Replace the rows with the corresponding rows from \(\mathbf{D}\).

8. The Chebyshev Helmholtz Solver¶

The 1D Helmholtz equation with Dirichlet BCs:

\[u'' - \alpha u = f, \quad u(-1) = a, \quad u(1) = b\]

In Chebyshev pseudospectral form, this becomes a (dense) linear system:

\[(\mathbf{D}^2 - \alpha \mathbf{I})\mathbf{u} = \mathbf{f}\]

with rows \(0\) and \(N\) replaced by the boundary conditions.

Rather than factorising this system on every call, spectraldiffx eliminates the two boundary unknowns (directly for Dirichlet, via the boundary rows of \(\mathbf{D}\) for Neumann) and diagonalises the remaining \((N-1)\times(N-1)\) interior operator once, \(\mathbf{E} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^{-1}\). \(\mathbf{E}\) depends only on the grid, so each solve is

\[\mathbf{u}_I = \mathbf{Q}\,(\boldsymbol{\Lambda} - \alpha)^{-1}\,\mathbf{Q}^{-1}\mathbf{r},\]

an \(O(N^2)\) operation for any \(\alpha\) — which may even be a traced JAX value. The eigenvalues are real and non-positive and \(\mathbf{Q}\) is well conditioned, so this is as accurate as a direct solve. For pure-Neumann Poisson (\(\alpha = 0\)) the constant eigenmode is dropped and the gauge \(u_{N/2} = 0\) is imposed.

For \(\alpha = 0\) this reduces to the Chebyshev Poisson solver.

9. Extension to 2D via Tensor Products¶

For 2D problems on \([-1,1]^2\), the Chebyshev basis extends naturally via tensor products:

\[u(x, y) = \sum_{m=0}^{M} \sum_{n=0}^{N} \hat{u}_{mn} T_m(x) T_n(y)\]

The 2D differentiation matrix for \(\partial/\partial x\) is:

\[\mathbf{D}_x = \mathbf{D} \otimes \mathbf{I}_N\]

and for \(\partial/\partial y\):

\[\mathbf{D}_y = \mathbf{I}_M \otimes \mathbf{D}\]

The 2D Laplacian is:

\[\nabla^2 = \mathbf{D}_x^2 + \mathbf{D}_y^2\]

For separable problems (e.g., rectangular domains with separable BCs), the 2D Helmholtz equation reduces to a sequence of 1D solves via the matrix diagonalisation method, reducing cost from \(O(N^4)\) to \(O(N^3)\). ChebyshevHelmholtzSolver2D uses exactly this (Haidvogel & Zang, 1979): the interior \(\mathbf{D}_x^2\) and \(\mathbf{D}_y^2\) blocks are diagonalised once, and each Dirichlet solve is four small matrix products plus a pointwise division by \(\lambda^y_j + \lambda^x_i - \alpha\).

References¶

Trefethen, L.N. (2000). Spectral Methods in MATLAB. SIAM.
Boyd, J.P. (2001). Chebyshev and Fourier Spectral Methods. Dover.
Canuto, C., et al. (2006). Spectral Methods: Fundamentals in Single Domains. Springer.
Haidvogel, D.B. & Zang, T. (1979). The accurate solution of Poisson's equation by expansion in Chebyshev polynomials. J. Comput. Phys. 30, 167–180.