Chebyshev and Multi-Domain Spectral Methods

Chebyshev methods approximate a smooth function globally rather than through local stencils. For analytic data they can gain many digits with modest polynomial order, but this rapid convergence is conditional: a corner, unresolved interface, poor coordinate map, or singular endpoint converts exponential decay into algebraic decay or worse. Multi-domain methods restore regularity by placing boundaries where the mathematical behavior changes and by enforcing the correct matching conditions.

Required background. Ordinary Differential Equations supplies boundary-value problems, interpolation, and eigenvalue equations.

Helpful background. Formulations of General Relativity supplies relativistic equations solved on radial and angular domains, while Partial Differential Equations and Green Functions clarifies regularity, boundary data, and separated operators.

Chebyshev collocation and mappings

On $-1\le x\le1$, choose Lobatto points

$$x_j=\cos\frac{j\pi}{N}, \qquad j=0,\ldots,N,$$

and represent $u_N$ by the degree-$N$ interpolating polynomial. Differentiation at the nodes becomes matrix multiplication, $u_N'=D_Nu_N$. The off-diagonal entries and a cancellation-preserving diagonal prescription are

$$(D_N)_{ij}=\frac{c_i}{c_j}\frac{(-1)^{i+j}}{x_i-x_j} \quad(i\ne j), \qquad (D_N)_{ii}=-\sum_{j\ne i}(D_N)_{ij},$$

where $c_0=c_N=2$ and the other $c_j$ equal one. The node ordering here runs from $+1$ to $-1$; reversing only the node array without transforming the matrix reverses or scrambles the derivative. Before solving an ODE, check $D_N\mathbf1=0$ and $D_N\mathbf x=\mathbf1$. For $N\ge2$, also check $D_N^2\mathbf x^2=2\mathbf1$. These polynomial identities catch implementation mistakes at negligible cost.

A physical interval $a\le r\le b$ follows from an affine map; a semi-infinite interval can use a rational map such as

$$r=L\frac{1+x}{1-x},$$

with endpoint regularity imposed on appropriately rescaled variables. The node $x=1$ represents infinity; do not evaluate the divergent map there and feed the result into an ordinary interior row. Derive a limiting endpoint equation for the rescaled variable. Polynomial approximation and collocation are developed in Grandclément and Novak 2009, §§ 2.1 and 2.4–2.5; rational compactification is treated separately in § 3.1.2, especially Eq. 89.

The scale $L$ is numerical input. If it places nearly all nodes far from the feature of interest, increasing $N$ can produce an ill-conditioned matrix before resolving the solution. Vary $L$ as part of the convergence study rather than tuning it once and omitting it from the record.

A regularity benchmark

Consider

$$u(x)=x\lvert x\rvert =\begin{cases}-x^2,&x\lt 0,\\x^2,&x\ge0.\end{cases}$$

It is continuously differentiable but its second derivative jumps at the origin. A single Chebyshev expansion on $[-1,1]$ therefore converges algebraically. Solve the manufactured boundary problem $u''=2\,\operatorname{sgn}x$ away from $x=0$, with $u(-1)=-1$ and $u(1)=1$, at $N=8,16,32,64$. Measure the maximum solution error on an independent uniform 4001-point grid, the differential residual outside a one-node neighborhood of the interface, the boundary residual, and the decay of Chebyshev coefficients.

Now split the interval into $[-1,0]$ and $[0,1]$. Each restriction is a quadratic polynomial, so degree two is sufficient in exact arithmetic, provided the interface conditions

$$u(0^-)=u(0^+), \qquad u'(0^-)=u'(0^+)$$

are enforced. The experiment cleanly demonstrates that spectral convergence belongs to a function together with its domain decomposition. Multi-domain collocation and matching are discussed in Grandclément and Novak 2009, §§ 2.6.1–2.6.4, whereas noncompact intervals use the distinct compactification construction in § 3.1.2.

Interfaces, tau rows, and conditioning

Boundary and interface equations replace selected collocation rows or enter through a tau formulation. Count conditions before solving: two second-order subdomains need four independent conditions. The two outer values plus continuity of $u$ and $u'$ supply exactly four. Omitting derivative continuity leaves a free homogeneous mode; imposing an extra condition overdetermines the system.

Map derivatives consistently. If $J=dr/dx$, then

$$\frac{d}{dr}=J^{-1}\frac{d}{dx}, \qquad \frac{d^2u}{dr^2}=J^{-2}u_{xx}-J^{-3}J_xu_x.$$

The second term disappears only for an affine map. The physical derivative, not the computational derivative, is matched at an interface. For a flux-form equation $(a(r)u_r)_r=s$, integrating through a source-free interface gives continuity of $a u_r$; if $a$ jumps, $u_r$ generally does not remain continuous. Derive the matching conditions from the equation rather than copying the constant-coefficient case.

Monitor both the matrix condition estimate and the continuous residual. Spectral differentiation matrices become ill-conditioned as $N$ grows even when the approximation improves. Rescaling rows, using integration formulations, and increasing arithmetic precision can separate representation conditioning from a genuine lack of convergence. Trefethen’s practical construction and boundary-row replacements are given in Trefethen 2000, Chapters 6–9, pp. 51–100.

Coefficient plots are diagnostic only when their normalization and parity are understood. For a smooth even or odd solution, one parity family may vanish exactly; do not mistake those zeros for the decay envelope. Plot the magnitude of the nonzero coefficients, the last resolved coefficient, and the dense-grid residual together. If coefficients fall geometrically and then level off at the arithmetic floor while the residual and solution error also saturate, precision is limiting. If the tail rises or oscillates before reaching that floor, the order, map, or domain split is inadequate.

Reproducibility and failure criteria

Pin the numerical library, transform normalization, node ordering, differentiation matrices, domain maps, interface enforcement, linear solver, arithmetic precision, and residual norms. A clean rerun evaluates the single-domain sequence above and two-domain degrees $N=2,4,8$ per subdomain, exports coefficients and condition estimates, samples the independent residual grid, and repeats the finest cases with 100-bit arithmetic. For the separate compactification check, use $L=1/2,1,2$ rather than a tuned single value. Exact reproduction should show algebraic single-domain behavior and near-roundoff multi-domain error once each quadratic is resolved.

The validity regime assumes smooth or analytic data inside each subdomain and correctly represented endpoint behavior. Failure modes include coefficient decay that stalls, endpoint contamination, a large off-grid residual despite nodal agreement, an omitted interface condition, and false convergence caused by precision saturation. Reject any claim based only on successive solutions agreeing: two under-resolved interpolants can agree while both miss the same localized feature.

Check your understanding

The two-domain calculation enforces both outer boundary values and continuity of $u$ at $x=0$, but not continuity of $u'$. Its nodal differential residual is zero. Why is the solution not determined?

Solution

Each subdomain solves a second-order equation, so four independent conditions are required. The two outer values and one continuity equation provide only three. A free homogeneous mode remains and can produce a derivative jump at the interface. Enforce the physically derived derivative or flux condition, then verify the interface residual and dense-grid convergence.

References

• Grandclément, P., and Novak, J. (2009), “Spectral Methods for Numerical Relativity,” Living Reviews in Relativity 12, 1. DOI.
• Trefethen, L. N. (2000), Spectral Methods in MATLAB, SIAM. DOI.