When integration sparsification fails: Banded Galerkin discretizations for Hermite functions, rational Chebyshev functions and sinh-mapped Fourier functions on an infinite domain, and Chebyshev methods for solutions with C${}^{\infty }$ endpoint singularities

Abstract

Chebyshev polynomial spectral methods are very accurate, but are plagued by the cost and ill-conditioning of dense discretization matrices. Modified schemes, collectively known as “integration sparsification”, have mollified these problems by discretizing the highest derivative as a diagonal matrix. Here, we examine five case studies where the highest derivative diagonalization fails. Nevertheless, we show that Galerkin discretizations do yield banded matrices that retain most of the advantages of “integration sparsification”. Symbolic computer algebra greatly extends the reach of spectral methods. When spectral methods are implemented using exact rational arithmetic, as is possible for small truncation $N$ in Maple, Mathematica and their ilk, roundoff error is irrelevant, and sparsification failure is not worrisome. When the discretization contains a parameter $L$, symbolic algebra spectral methods return, as answer to an eigenproblem, not discrete numbers but rather a plane algebraic curve defined as the zero set of a bivariate polynomial $P(\lambda ,L)$; the optimal approximations to the eigenvalues ${\lambda }_j$ are in the middle of the straight portions of the zero contours of $P(\lambda ;L)$ where the isolines are parallel to the $L$ axis.

Introduction

Single domain spectral methods generate dense discretization matrices. However, Charles Clenshaw showed in 1957 that for some special differential equations, Chebyshev–Galerkindiscretizations are sparse. Berioz, Hagstrom, Lau and Price [2] coined the neologism “integration sparsification” for the numerical strategies that discretize the highest derivative as a diagonal matrix or a matrix of very low bandwidth. (A matrix is “banded” with bandwidth $m$ if its elements have the property that $G_{j,k}=0$ for all $|j-k|>m$.)

Integration sparsification is reviewed at length in [26]. There are multiple strategies that seemingly have very different rationales and mechanics, but the end result is the same or almost the same for all. One motive for this recent revival of “integration sparsification” is that standard pseudospectral methods are somewhat ill-conditioned for large truncation $N$ and/or high order differential equation. Integration sparsification greatly reduces the condition number of the discretization matrix and also facilitates fast direct methods for solving the linear algebra problems.

In this article, we analyze five case studies where integration sparsification fails. The connections between them will at times seem rather loose. Much of the text is not about great ideas, but rather about the craft knowledge, the little tricks and skills that distinguish a master brewer making consistent, reliable beers, ales and largers from an amateur zymurgist making beer in a bathtub. Collectively, though, these examples support a number of important themes:

• Highest derivative diagonalization is often impossible.

• Nevertheless, Galerkin discretizations do yield banded matrices that retain most of the advantages of “integration sparsification”.

• Symbolic computer algebra greatly extends the reach of spectral methods. When spectral methods are implemented using exact rational arithmetic, as is possible for small truncation $N$ in Maple, Mathematica and their ilk, roundoff error is irrelevant, and sparsification failure is not worrisome.

• The most sparsified algorithm is not necessarily the best; the rate of convergence cannot be ignored. In particular sinh-mapped domain truncation is usually superior to unmapped Fourier truncation even though only the latter yields diagonal derivative discretizations.

• Symbolic spectral methods return an algebraic plane curve $P(\lambda ;L)=0$ as the numerical solution to an eigenvalue problem where $\lambda$ is the eigenvalue and $L$ is a numerical parameter; $L$ is optimum where the contours are straight lines parallel to the $L$-axis.

• Integration sparsification also greatly facilitates the preconditioned matrix iteration schemes proposed more than 40 years ago by L. M. Delves and his collaborators  [17], [18], [19]; this is a major theme of our articles [24], [25], [26], [27] and so will receive only a limited discussion.

The Delves preconditioner operates in “modal space” (that is on the Galerkin matrix) rather than on the nodal basis and collocation matrix. Instead of solving a matrix problem with the full, dense spectral discretization, Delves proposed iterating with a matrix preconditioner that is very sparse. Galerkin approximations are usually not “diagonally dominant” in the classical sense. However, for a Fourier Galerkin matrix or a “sparsified” Chebyshev discretization, the diagonal element of the $n$th row of the matrix does become more and more important relative to the other elements in the same row and column as $n\to \infty$. In this sense, the Galerkin matrix is “asymptotically diagonal”. In our other works [24], [25], [26], [27], we focus on the merits of preconditioned iterations that exploit the “asymptotic diagonality” for “sparsified” spectral methods.

The iterations are needed because genericallly spectral matrices are dense and direct, non-iterative methods have costs that scale as the cube of the number of degrees of freedom $N$. In contrast, preconditioned iterations reduce the cost by orders of magnitude, but only if a cheap, effective preconditioner is available. “Sparsified” derivative matrices are valuable because they are excellent preconditioners. They drastically reduce condition number because they discretize the highest order derivative which is the primary source of the ill-conditioning in spectral discretizations; they are cheap to apply as preconditioners because of their sparsity.

To explain why sparsification sometimes fails, we need to explain what it is. In the next section, we explain the Petrov–Galerkin method in which the basis functions are Chebyshev polynomials but the “test functions” are ultraspherical (Gegenbauer) polynomials.

There is a spectrum of options for populating the matrix elements. Matrix–matrix multiplication is easy-to-program but computationally expensive. Other options trade chirugery¹ for fewer floating point operations as explained in Section 3.

We then show that for problems on an infinite domain, it is usually not possible to diagonalize the discretization of the highest derivative in a differential equation. For both Hermite functions and rational Chebyshev functions, however, we show that it is easy to generate a banded discretization.

Next, we compare Fourier domain truncation for an unbounded interval applied both with and without the “sinh-map”. Ironically, the unmapped method which trivially generates diagonal derivative matrices not just for the highest derivative but for all derivatives, is much less efficient than the sinh-map variant, whose derivatives are not sparsified.

We then turn to problems on a finite interval. We show that for differential equations that are singular at an endpoint, the standard sparsification procedures cannot be applied. Nevertheless, for solutions which are not analytic but are infinitely differentiable even at the singularity, i.e., belong to the $C^{\infty }$ function class, it is still possible to obtain an exponential rate of convergence from a banded discretization of low bandwidth.