Original articlesWhen 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 endpoint singularities
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 if its elements have the property that for all .)
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 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:
- 1.
Highest derivative diagonalization is often impossible.
- 2.
Nevertheless, Galerkin discretizations do yield banded matrices that retain most of the advantages of “integration sparsification”.
- 3.
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 in Maple, Mathematica and their ilk, roundoff error is irrelevant, and sparsification failure is not worrisome.
- 4.
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.
- 5.
Symbolic spectral methods return an algebraic plane curve as the numerical solution to an eigenvalue problem where is the eigenvalue and is a numerical parameter; is optimum where the contours are straight lines parallel to the -axis.
- 6.
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 th row of the matrix does become more and more important relative to the other elements in the same row and column as . 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 . 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 chirugery1 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 function class, it is still possible to obtain an exponential rate of convergence from a banded discretization of low bandwidth.
Section snippets
Petrov–Galerkin discretizations
We assume most readers have at least a nodding acquaintance with Galerkin discretizations and Chebyshev polynomial methods, so this brief review is mostly to establish notation and terminology. For the novice, good introductions are [11], [40].
Choose a set of basis functions , a set of test functions and the weighted inner product where is the weight function. Write a linear problem as where is a linear differential or integral operator. Then
Generating the Galerkin matrix
Pseudospectral algorithms, which evaluate rather than integrate, became dominant because of their simplicity. However, once the pseudo-spectral matrix has been computed, generation of the Galerkin matrix is easy.
Let denote the pseudospectral matrix with elements where is the operator of the differential equation and let be a vector of samples of at the quadrature points.
Compute the quadrature weights and evaluate test function at each point
Hermite functions
The Hermite functions arise as the solutions to a wide variety of physics problems. They are the eigenfunctions of the quantum harmonic oscillator; they are factors in the separation-of-variables solution to PDEs in parabolic cylinder coordinates; they are the latitudinal structure of planetary scale equatorial waves in the ocean and atmosphere, very important in climate. In addition, they are a good general basis for solving problems on an infinite domain, provided that the solutions decay
Fourier domain truncation and the sinh-mapping
Another alternative for an unbounded interval is “Fourier domain truncation”, which is a linear stretching transformation from the computational coordinate to the physical coordinate accompanied by total disregard for everything outside the truncated domain in . The chain rule generates the derivative transformations: An essential prerequisite to obtaining an exponential convergence rate without additional tricks is that must decay
Computer algebra and spectral methods
Symbolic manipulation systems like Maple, Mathematica, REDUCE and the Symbolic Toolbox in Matlab are powerful partners with spectral methods [1], [9], [13], [34]. Computer algebra systems can compute Galerkin matrices for many problems exactly. For example, the eigenvalues of the quantum quartic oscillator are, for a basis truncated at Hermite functions, the roots of the determinant of the matrix:
solutions to ODEs and the failure of ultraspherical diagonalization
Cornelius Lanczos, who was the pioneer of solving differential equations using Chebyshev polynomials, triumphed, even in his first article (1938) [30], over ordinary differential equations with endpoint singularities and solutions that were infinitely differentiable (that is, belong to the function class ) but not holomorphic at an endpoint. He showed that the divergent asymptotic reciprocal power series for the exponential integral could be replaced easily by an exponentially convergent
Summary and open problems
We have shown through a mixture of analysis and examples that diagonalization of the highest derivative is not always possible for spectral methods. In each case, however, we still obtain a sparse, banded discretization for differentiation matrices. The discretization matrix of a differential equation, which typically involves transcendental functions multiplying derivatives, are usually full, dense matrices. Direct solution of a dense matrix is proportional to the cube of matrix size .
Acknowledgments
Support was from NSF grant OCE 1059703, a China Scholarship Council grant 201306280061 and NSFC grant key program (No. 51236006) of China.
References (43)
- et al.
Investigating a hybrid perturbation-galerkin technique using computer algebra
J. Symbolic Comput.
(1991) The asymptotic coefficients of Hermite series
J. Comput. Phys.
(1984)Spectral methods using rational basis functions on an infinite interval
J. Comput. Phys.
(1987)Multipole expansions and pseudospectral cardinal functions: A new generalization of the Fast Fourier Transform
J. Comput. Phys.
(1992)Chebyshev and Legendre spectral methods in algebraic manipulation languages
J. Symbolic Comput.
(1993)The rate of convergence of Fourier coefficients for entire functions of infinite order with application to the Weideman-Cloot sinh-mapping for pseudospectral computations on an infinite interval
J. Comput. Phys.
(1994)- et al.
An adaptive algorithm for spectral computations on unbounded domains
J. Comput. Phys.
(1992) - et al.
High-order shifted Gegenbauer integral pseudo-spectral method for solving differential equations of Lane-Emden type
Appl. Numer. Math.
(2018) - et al.
Modal preconditioning of Galerkin spectral methods: Dual bookkeeping for the Delves-Freeman iteration
J. Comput. Phys.
(2015) - et al.
Bandwidth Truncation for Chebyshev polynomial and Ultraspherical/Chebyshev Galerkin discretizations of differential equations: Restrictions and two improvements
J. Comput. Appl. Math.
(2016)