Construction of Laguerre pseudospectral differentiation matrices

1. Problem Statement

Laguerre pseudospectral methods are essential for solving differential equations on semi-infinite domains $[0, \infty)$. However, the construction of Laguerre pseudospectral differentiation matrices (D) suffers from severe numerical instability as the number of collocation points ($n$) increases.

• The Core Issue: Classical formulations involve intermediate quantities (such as Laguerre polynomials $L_n(x)$ and their derivatives, or barycentric weights $v_j$) that grow or decay exponentially. In double-precision floating-point arithmetic, these intermediate terms exceed the representable range (overflow or underflow) even when the final matrix entries are well-scaled.
• Consequence: Standard implementations, such as the widely used DMSUITE package (Weideman, 2003), experience catastrophic numerical breakdown at relatively modest values of $n$ (approximately $n \approx 125$). This limits the applicability of Laguerre methods for high-accuracy or high-resolution problems.
• Specific Instabilities:
  • Overflow: $L_n(x)$ grows rapidly for large $x$.
  • Underflow: The associated Laguerre functions $\tilde{L}_n(x) = e^{-x/2}L_n(x)$ require evaluating $e^{-x/2}$, which underflows for large $x$.
  • Catastrophic Cancellation: Standard formulas for diagonal entries often rely on the "negative sum trick" (summing off-diagonal entries), which amplifies round-off errors in weighted interpolation settings.

2. Methodology

The authors propose a stable, "all-in-one" construction procedure that reformulates the matrix entries and computes them simultaneously to avoid intermediate overflow/underflow.

A. Reformulation of Off-Diagonal Entries

Instead of computing the unscaled coefficients $c_j = e^{-x_j/2}/v_j$ separately (which leads to overflow in $1/v_j$ and underflow in $e^{-x_j/2}$), the authors introduce scaled coefficients $\tilde{c}_j$:
$$\tilde{c}_j = a(x_j) \tilde{L}_n^{(\alpha)\prime}(x_j)$$
where $\tilde{L}_n^{(\alpha)\prime}$ is the derivative of the Laguerre function, and $a(x)$ depends on the node distribution (1 for standard, $x$ for augmented/Radau).

• Key Insight: The off-diagonal entries $D_{k,j}$ depend only on the ratio $\tilde{c}_k / \tilde{c}_j$. By computing the numerator and denominator simultaneously, the exponential factors cancel out naturally without ever being explicitly evaluated at extreme values.

B. The Glaser–Liu–Rokhlin (GLR) Algorithm

To compute the required roots (collocation nodes) and the derivatives of the Laguerre functions simultaneously, the authors utilize the Glaser–Liu–Rokhlin algorithm (Glaser et al., 2007).

• Mechanism: This algorithm solves the second-order ODE satisfied by Laguerre functions using a predictor-corrector approach. It advances from a known root $x_i$ to $x_{i+1}$ using local Taylor expansions derived from the ODE.
• Stability: Because the algorithm works directly with the Laguerre functions (which include the $e^{-x/2}$ weight) and uses local expansions, it avoids the explicit evaluation of large polynomials or small exponentials.
• Modification for Small $x$: Near the origin ($x=0$), the ODE has a singularity causing dense root clustering. The authors adapt the algorithm by using a modified recurrence relation (Huang & Yu, 2024) for the first few roots to ensure accuracy before switching to the Taylor expansion method.

C. Closed-Form Expression for Diagonal Entries

The authors reject the "negative sum trick" for diagonal entries in favor of a direct closed-form formula derived from the properties of Laguerre polynomial zeros (Stieltjes relations).

• Standard Nodes: $D_{k,k} = -\frac{1+\alpha}{2x_k}$.
• Augmented/Radau Nodes: Specific formulas are derived that account for the inclusion of $x=0$.
• Advantage: This approach avoids the accumulation of round-off errors inherent in summing many off-diagonal terms and allows for the computation of individual diagonal entries without needing all other nodes (useful for truncated spectral methods).

3. Key Contributions

1. Stable Reformulation: A new mathematical formulation for differentiation matrix entries that eliminates catastrophic cancellation and overflow/underflow by utilizing scaled coefficients and simultaneous computation.
2. Unified Algorithm: Integration of the Glaser–Liu–Rokhlin algorithm to generate both collocation nodes and the necessary derivative values in a single, stable pass.
3. Direct Diagonal Formulas: Derivation and validation of explicit closed-form expressions for diagonal entries that outperform the traditional negative sum trick in terms of accuracy and independence from other nodes.
4. Comprehensive Implementation: A complete procedure covering standard, augmented, and Radau node distributions for both first and second-order differentiation matrices.

4. Results

Numerical experiments were conducted comparing the proposed method against the standard DMSUITE implementation (using 24-digit arithmetic as a reference).

• Robustness: The proposed method remains stable and accurate for $n$ up to 500 and beyond, whereas DMSUITE fails (returns Inf or NaN) at $n \approx 125$.
• Accuracy:
  • Off-diagonal entries: The proposed method achieves machine precision ($\sim 10^{-16}$) for large $n$, while DMSUITE errors grow rapidly.
  • Diagonal entries: The direct formula yields errors comparable to the root-finding accuracy, significantly outperforming the negative sum trick used in DMSUITE.
  • Second-order matrices: The proposed method maintains high accuracy for second derivatives, where DMSUITE breaks down due to instability in the underlying first-order matrix.
• Application Examples:
  • Schrödinger Equation: The method successfully computed eigenvalues for the Woods–Saxon potential with exponential convergence, resolving high-order eigenstates that were inaccessible to DMSUITE.
  • Boundary Value Problems: A model problem on the half-line was solved with near-machine precision ($10^{-15}$) using $n \approx 230$ points, a regime where standard methods fail.

5. Significance

This work resolves a long-standing practical limitation in spectral methods on semi-infinite domains.

• Extending Applicability: It enables the use of Laguerre pseudospectral methods for problems requiring high resolution (large $n$), which was previously impossible due to numerical breakdown.
• Reliability: By providing a stable, "all-in-one" construction, it removes the need for ad-hoc scaling or complex stabilization techniques that often degrade accuracy.
• Foundation for Future Work: The stable construction of differentiation matrices and nodes paves the way for advanced techniques such as truncated spectral collocation, adaptive scaling, and efficient interpolation on unbounded domains.

In summary, Nel and Hale provide a mathematically rigorous and numerically robust solution to the construction of Laguerre differentiation matrices, effectively modernizing the toolkit for solving differential equations on $[0, \infty)$.