---
title: Chebyshev Pseudospectral Method Overview
url: https://www.emergentmind.com/topics/chebyshev-pseudospectral-method-cpm
type: topic
---

# Chebyshev Pseudospectral Method Overview

Chebyshev Pseudospectral Method (CPM) denotes a family of spectral discretizations in which functions are approximated by Chebyshev polynomials and differential operators are enforced either at selected collocation nodes or in coefficient space. In the cited literature, CPM encompasses Chebyshev–Gauss and Chebyshev–Gauss–Lobatto nodal schemes, shifted-Chebyshev operational-matrix constructions on finite intervals, time-space collocation, and the Chebyshev-coefficient “ultraspherical” method of Sheehan Olver and Alex Townsend for variable-coefficient ODEs [1202.1347]. Across these variants, the common structure is the replacement of continuous derivatives, products, and boundary conditions by algebraic operators with spectral accuracy for sufficiently smooth data, together with problem-specific devices for conditioning, adaptivity, and the treatment of singular sources, nonhomogeneous data, fractional operators, and optimal-control constraints.

## 1. Basis functions, grids, and interval mappings

The foundational object is the Chebyshev polynomial of the first kind,
\[
T_n(x)=\cos\!\bigl(n\arccos x\bigr),\qquad x\in[-1,1],
\]
with the three-term recurrence
\[
T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).
\]
Several formulations also use the second-kind family \(U_n\), especially when derivatives or endpoint conditions are built directly into the basis [1202.1347; 2505.19454].

Two nodal sets dominate the literature. The Chebyshev–Gauss–Lobatto nodes are
\[
x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,
\]
while the Chebyshev–Gauss nodes are
\[
x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.
\]
The cited papers also use the sign-reversed Lobatto ordering
\[
x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,
\]
which is the same grid with opposite orientation [1802.03405]. On \([0,1]\), the affine map
\[
x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]
\]
produces Chebyshev–Gauss–Lobatto nodes adapted to shifted polynomials and diffusion problems on the unit square [2509.19338]. On a general interval \([a,b]\), one uses
\[
x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,
\]
and for moving time windows \([t_k,t_{k+1}]\) the shifted basis is
\[
\phi_n(t)=T_n\!\Bigl(\frac{2(t-t_k)}{\Delta t_k}-1\Bigr),\qquad \Delta t_k=t_{k+1}-t_k,
\]
with collocation instants obtained by mapping canonical Chebyshev nodes into physical time [2505.07234].

Two approximation viewpoints recur. In coefficient-space methods, one writes
\[
u(x)\approx \sum_{k=0}^n u_k\,T_k(x),
\]
so the unknown is the Chebyshev coefficient vector. In nodal collocation, one instead interpolates by Lagrange polynomials at Chebyshev nodes; the time-space Burgers–Fisher scheme, for example, uses the one-dimensional Chebyshev Lagrange basis
\[
\ell_k(z)=\prod_{\substack{m=0\\m\neq k}}^N\frac{z-z_m}{z_k-z_m}
\]
at Chebyshev–Gauss–Lobatto points in both \(x\) and \(t\) [2306.09988]. The shifted-Chebyshev fractional formulation likewise expands
\[
\widetilde T_n(t)=T_n(2t-1)=\cos\bigl(n\arccos(2t-1)\bigr)
\]
on \([0,1]\) [2511.11082].

## 2. Differentiation matrices, multiplication operators, and assembled systems

In nodal collocation, differentiation is realized by dense differentiation matrices. On a Chebyshev–Lobatto grid with \(c_0=c_N=2\) and \(c_j=1\) for \(1\le j\le N-1\), the first-derivative matrix has entries
\[
D^{(1)}_{jk}=\frac{c_j(-1)^{j+k}}{c_k(x_j-x_k)},\qquad j\neq k,
\]
with diagonal entries defined by the row-sum identity
\[
D^{(1)}_{jj}=-\sum_{k\neq j} D^{(1)}_{jk}.
\]
Second derivatives are formed either as \(D^{(2)}=D^{(1)}D^{(1)}\) or from the explicit Lobatto formula [1802.03405]. In time-space collocation, these one-dimensional operators are lifted to tensor form through Kronecker products:
\[
\frac{\partial}{\partial t}\longleftrightarrow I_x\otimes D_t,\qquad
\frac{\partial}{\partial x}\longleftrightarrow D_x\otimes I_t,\qquad
\frac{\partial^2}{\partial x^2}\longleftrightarrow D_x^2\otimes I_t
\]
[2306.09988].

The coefficient-space ultraspherical method uses a different algebra. If
\[
a(x)=\sum_{j=0}^{m-1} a_j\,T_j(x),\qquad
u(x)=\sum_{k=0}^n u_k\,T_k(x),
\]
then multiplication by \(a(x)\) becomes a Toeplitz \(+\) Hankel operator \(M[a]\) of bandwidth \(m\). Higher derivatives are represented by diagonal operators \(D_\lambda\) in the ultraspherical basis \(C^{(\lambda)}\), and conversion operators \(S_\lambda\) map \(C^{(\lambda)}\)-coefficients to \(C^{(\lambda+1)}\)-coefficients. The \(\lambda\)th term of a variable-coefficient ODE therefore acts as
\[
S_{N-1}\cdots S_\lambda\,D_\lambda\,M[a^\lambda]
\]
on the Chebyshev-coefficient vector [1202.1347].

A central structural consequence is “almost bandedness.” After imposing \(K\) boundary conditions by row replacement and permuting them to the top, the global matrix contains \(K\) dense boundary rows, a fixed number of diagonals from differentiation and conversion, and a multiplication band of width \(m\). Because \(m\) and the differential order \(N\) are fixed by the smoothness of the coefficients, the discretized system is almost banded rather than fully dense [1202.1347].

The same operator logic extends to higher dimensions and mixed bases. In the transient anisotropic diffusion formulation on \(\Omega=(0,1)\times(0,1)\), the full Chebyshev differentiation matrix \(D\) is sliced into interior-only blocks \(D_1\) and \(D_2\), and anisotropy enters through diagonal matrices of sampled conductivities together with boundary corrections derived from Robin-type conditions [2509.19338]. In superconducting strip problems on \(x\in\mathbb R\), \(y\in[-1,1]\), the method combines Chebyshev differentiation in the finite \(y\)-direction with Hermite-function differentiation in the unbounded \(x\)-direction, producing
\[
D_y=RCP^{-1},\qquad D_x=\frac1L\,Q\,B\,Q^{-1}
\]
for the two coordinate directions [2108.10654].

## 3. Conditioning, preconditioning, and solver architecture

The best-developed conditioning theory in the cited set is the Olver–Townsend ultraspherical scheme for linear ODEs with variable coefficients. The key diagonal preconditioner is
\[
P=\operatorname{diag}(1,\dots,1,1/N,1/(N+1),\dots,1/n),
\]
which balances the growth of ultraspherical differentiation and conversion so that
\[
PA=I+K
\]
with \(K\) compact in an appropriate \(\ell^2_\lambda\) norm [1202.1347]. For Dirichlet boundary conditions, one may take \(\lambda=0\), yielding a bounded \(2\)-norm condition number of \(PA\) uniformly in \(n\). In that sense, the method is well conditioned in the standard Euclidean norm for that class of boundary data.

The associated direct solver is an adaptive QR factorization applied to an “infinite” almost-banded operator. Givens rotations are introduced column by column, each filled-in row is maintained in \(O(m+K)\) storage, and truncation is determined by the transformed right-hand side. The stopping test is
\[
\|r_{j+1:}\|_2\le \varepsilon |r_1|,
\]
with a representative tolerance \(10^{-14}\). The operation count is \(O(nm^2)\), the memory requirement is \(O(nm)\), and the reported implementation can “efficiently and reliably solve for solutions that require as many as a million unknowns” [1202.1347].

Other CPM variants emphasize stability through different mechanisms. The time-space Burgers–Fisher scheme derives a discrete weighted-norm estimate,
\[
\|U(\cdot,t)\|_w\le \exp(Ct)\,\|U(\cdot,0)\|_w,
\]
after discrete integration by parts, thereby proving spectral stability of the fully coupled time-space collocation method [2306.09988]. The shifted-Chebyshev fractional construction identifies the dominant numerical hazard not in the discrete-cosine transform, whose conditioning is \(O(1)\), but in the coefficient matrix \(C\) used to generate operational matrices. Because its entries grow like \(O(4^N)\) and alternate in sign, the method uses symbolic or variable-precision arithmetic with digits chosen approximately linearly in \(N\) to preserve stability up to \(N=1000\) in the reported tests [2511.11082].

## 4. Boundary conditions, singular sources, and PDE generalizations

A persistent theme in CPM research is that singular or nonhomogeneous data are typically not inserted naively into a single global discretization. The “Particle-without-Particle” method is a direct illustration. For a one-dimensional linear PDE with a source
\[
L[u](x)=\sum_{\ell=0}^m a_\ell(x)\,u^{(\ell)}(x)=\sum_{k=0}^K \alpha_k\,\delta^{(k)}(x-x_p),
\]
the domain is split into \(\mathcal D^-=[a,x_p]\) and \(\mathcal D^+=[x_p,b]\), each with its own Chebyshev–Lobatto grid. One solves the homogeneous PDE \(L[u^\pm]=0\) in each subdomain and replaces the singular source by jump conditions at the interface. For the model problem
\[
u''+u=a\,\delta(x)+b\,\delta'(x),
\]
the jumps are
\[
[u]=b,\qquad [u']=a
\]
[1802.03405]. The paper proves that the method applies to any linear PDE whose source is any linear combination of delta distributions and derivatives thereof supported on a one-dimensional subspace of the problem domain.

For nonlinear initial-boundary-value problems, nonhomogeneous data are often homogenized before collocation. In the time-space Chebyshev pseudospectral method for the generalized Burgers–Fisher equation,
\[
U_t + A\,(U^{\delta+1})_x - B\,U_{xx} - C\,(U-U^{\delta+1})=0,
\]
a lifting function \(\Omega(x,t)\) interpolates the initial and boundary data, and the unknown is rewritten as
\[
U(x,t)=V(x,t)+\Omega(x,t)
\]
so that \(V\) satisfies homogeneous initial and boundary conditions [2306.09988]. Interior collocation then yields a nonlinear algebraic system of dimension \((N-1)\times N\), solved by Newton–Raphson.

In transient anisotropic diffusion, the spatial operator is built from a time-dependent rotated diffusion tensor
\[
K(t,x,y)=\Phi(t)\Lambda(x,y)\Phi^T(t),\qquad
\Lambda=\operatorname{diag}[k_{11}(x,y),k_{22}(x,y)],
\]
with \(\Phi(t)\) determined by a rotation angle \(\theta(t)\). Robin-type conditions on all four sides of the square are reformulated to express normal derivatives in terms of \(U\) and boundary data, after which lexicographic ordering yields the semidiscrete system
\[
U'(t)=M(t)\,U(t)+Sg(t),\qquad U(0)=U_0.
\]
Time advancement is performed by Crank–Nicolson, requiring a solve with \(I-(\Delta t/2)M\) at each step when \(M\) varies in time [2509.19338].

The Hermite–Chebyshev superconducting-strip formulation shows a further geometric generalization. On \(x\in\mathbb R\), \(y\in[-1,1]\), the unknown is expanded in a tensor basis of Hermite functions and Chebyshev polynomials, while a stream function \(g\) is expanded in second-kind Chebyshev polynomials \(U_m(y)\) so that the Dirichlet condition \(g(\pm 1)=0\) is built into the basis. At each time step, the method computes current density from spectral derivatives, evaluates a nonlinear current-voltage law pointwise, solves \(N+1\) dense systems of size \(M+1\), reconstructs the field, and advances with an ODE solver such as MATLAB’s `ode15s` [2108.10654].

## 5. Fractional operators, optimization, control, and online identification

The CPM framework has also been extended well beyond classical integer-order boundary-value problems. For Caputo-type advection-diffusion equations, shifted Chebyshev polynomials on \([0,1]\) are used to construct an operational matrix \(D^{(\alpha)}\) for the Caputo derivative and an operational matrix \(I^{(\alpha)}\) for the Riemann–Liouville integral. If
\[
f(t)\approx \sum_{k=0}^N \hat f_k\,\widetilde T_k(t/T),
\]
then at the shifted Chebyshev–Gauss–Lobatto nodes \(t_j=\frac{T}{2}[1+\cos(j\pi/N)]\),
\[
f^{(\alpha)}(t_j)\approx \sum_{k=0}^N D^{(\alpha)}_{j+1,k+1}\,\hat f_k.
\]
For the one-dimensional PDE
\[
\partial_t^\alpha u=\kappa\,\partial_x^2u + v\,\partial_xu + f(x,t),
\]
the discretization leads to a Sylvester equation
\[
A\,\text{inner}+\text{inner}\,B=C,
\]
solved by MATLAB’s built-in `lyap`; in \(d\) spatial dimensions it becomes a tensor Sylvester equation solved by `sylvesterND` from Cuesta–de la Hoz ’24 [2511.11082].

In direct optimal control, the cited methodology transcribes an ODE-constrained problem into a nonlinear program by approximating only the highest-order derivative in a Chebyshev series and recovering lower-order states by successive integration. If \(x^{(m)}\) is approximated as
\[
x^{(m)}(t(\tau))\approx \sum_{j=0}^N c_j\,T_j(\tau),
\]
then sparse integration matrices \(B^{(1)},\dots,B^{(m)}\) map the coefficient vector \(c\) to lower derivatives and states. Dynamic constraints are collocated at Chebyshev nodes, endpoint constraints are enforced directly for CGL or by interpolation for CG, and the cost is integrated by Clenshaw–Curtis weights. The decision vector may take the form
\[
\chi=[c_0,\dots,c_N;\,u_0,\dots,u_N;\,t_f]
\]
[2505.19454].

A related but distinct line of work uses CPM for one-dimensional optimization. The exact line-search method of “Optimization via Chebyshev Polynomials” constructs a fourth-order Chebyshev interpolant on five CGL points, differentiates it to obtain a cubic model, classifies its roots analytically, and then refines the candidate minimizer by a Newton step using first- and second-order Chebyshev pseudospectral differentiation matrices. If the Hessian surrogate is non-positive or the iterate exits the interval, the method falls back to a golden-section or Brent step. The paper states that the resulting update converges quadratically when second-order information is used [1603.00863].

The online identification work extends the method to moving time windows. Over \(I_w=(t^{w-1},t^w]\), the system drift is approximated as
\[
\dot x(t)\approx \sum_{n=0}^{N_w} a_n\,\phi_n(t),
\]
with \(N_w+1\) Chebyshev sampling instants. Coefficients are computed by the regularized least-squares formula
\[
a^w=(A^\top A+R_0)^{-1}(R_0a_0+A^\top b),
\]
and continuity across windows is imposed by matching all derivatives up to order \(N_w\) at \(t=t^w\). The same piecewise-Chebyshev approximation is then embedded in an observer whose gain solves a Lyapunov equation, and the paper proves uniform boundedness of the parameter and state errors [2505.07234].

## 6. Convergence behavior, quantitative results, and methodological scope

The dominant accuracy claim across the cited literature is spectral or exponential convergence for sufficiently smooth or analytic solutions. In the Particle-without-Particle formulation, if each subdomain solution \(u^\pm\) is smooth, then Chebyshev collocation yields
\[
E(N)=\|u-u_N\|=O(e^{-\sigma N})
\]
for some \(\sigma>0\) determined by analyticity-strip width [1802.03405]. The time-space Burgers–Fisher study reports max-norm errors decaying from \(2.27\times 10^{-9}\) at \(N=4\) to \(2.22\times 10^{-16}\) at \(N=16\), which is the expected spectral decay until round-off [2306.09988]. The anisotropic diffusion abstract states that analytic solutions exhibit \(O(e^{-\alpha n})\) accuracy, whereas solutions with only \(m\) continuous derivatives recover \(O(n^{-m})\) algebraic decay [2509.19338]. The Hermite–Chebyshev strip method is spectrally convergent in \(y\) but only algebraically convergent in \(x\) when the solution decays with a power law at infinity [2108.10654].

Several papers report concrete performance data.

| Context | Reported outcome | Source |
|---|---|---|
| Variable-coefficient ODE \(u'+x^3u=100\sin(2\times10^4 x^2)\), \(u(-1)=0\) | adapts to \(n\approx 2.0\times 10^4\) for 16-digit accuracy; condition number of \(PA\) stays \(O(1)\); solve time \(\propto n\) | [1202.1347] |
| Airy BVP with \(\varepsilon=10^{-9}\) | solver finds \(n\approx 2.0\times 10^4\) in \(\approx 10\,\mathrm{s}\) (C++), with uniform error \(\lesssim 10^{-12}\) | [1202.1347] |
| Advection equation with moving \(\delta\) source | at \(N=160\), PwP + Chebyshev PSC gives \(L^2\) error \(\sim 10^{-4}\) versus \(\sim 10^{-3}\) for \(\delta\)-approximation + 6th-order finite difference | [1802.03405] |
| Elliptic Poisson problem with ring \(\delta\) source | at \(N=40\), PwP + CL-PSC gives \(L^\infty\) error \(\sim 10^{-6}\) versus \(\sim 10^{-3}\) for Tornberg–Kreiss \(\delta\)-approximation | [1802.03405] |
| Inhomogeneous superconducting strip | \((M,N)=(50,80)\) gives \(\|\Delta j\|_{L^\infty}\approx 1.6\%\); \((60,100)\) reduces error to \(\sim 0.6\%\) | [2108.10654] |
| Superconducting dynamo benchmark | voltage error \(\sim 0.9\%\) with \((M,N)=(40,60)\), CPU time \(\sim 3\)–\(5\) min/cycle | [2108.10654] |
| Rocket landing flip maneuver | with \(N=60\), Chebyshev CG/CGL/second-kind and Legendre LGL give final mass \(\approx 93\,980\) kg in \(\sim 13.58\) s; differences \(<0.05\) kg | [2505.19454] |
| One-dimensional optimization tests | CPSLSM achieved 10–15 correct digits versus MATLAB’s `fminbnd` typically 5–8 digits; iterations 2–15 versus 5–40 | [1603.00863] |

These results also clarify the scope of the method. The cited literature documents several counterexamples to the view that CPM is limited to a single dense differentiation matrix on one interval. The method appears as an almost-banded coefficient-space scheme for linear ODEs [1202.1347], a multi-domain collocation method with exact jump conditions for distributional sources [1802.03405], a fully coupled time-space collocation method [2306.09988], a shifted-Chebyshev operational-matrix construction for fractional derivatives [2511.11082], a mixed Hermite–Chebyshev method on \(\mathbb R\times[-1,1]\) [2108.10654], and an integral-collocation transcription for nonlinear programming in optimal control [2505.19454].

At the same time, the literature records clear limitations. If coefficient functions develop singularities near \([-1,1]\), the number \(m\) of modes required to represent them may grow with \(n\), and the direct ultraspherical solver remains stable but its cost becomes \(O(n^3)\) [1202.1347]. Highly oscillatory temporal data in the fractional setting require \(N\sim O(m)\) to resolve frequencies \(e^{im t}\) [2511.11082]. In infinite-domain strip problems, the Hermite direction loses spectral convergence when the underlying solution has only algebraic decay [2108.10654]. These caveats do not invalidate CPM; rather, they delimit the regimes in which its characteristic spectral efficiency is fully realized.

Source: https://www.emergentmind.com/topics/chebyshev-pseudospectral-method-cpm