Hybrid Chebyshev function bases for sparse spectral methods in parity-mixed PDEs on an infinite domain

TL;DR

This paper introduces a spectral method using hybrid Chebyshev function bases to efficiently solve differential equations on infinite domains, maintaining sparsity even with parity-mixed operators, enabling fast and accurate computations.

Contribution

The paper develops a novel hybrid Chebyshev rational function basis that preserves sparsity in spectral discretizations of parity-mixed PDEs on infinite domains.

Findings

Achieves $O(N \, \ln N)$ computational complexity.

Maintains spectral accuracy in unbounded geometries.

Enables fast implicit-explicit time-marching schemes.

Abstract

We present a numerical spectral method to solve systems of differential equations on an infinite interval $y\in (-\infty, \infty)$ in presence of linear differential operators of the form $Q(y) \left(\partial/\partial_y\right)^b$ (where $Q(y)$ is a rational fraction and $b$ a positive integer). Even when these operators are not parity-preserving, we demonstrate how a mixed expansion in interleaved Chebyshev rational functions $TB_n(y)$ and $SB_n(y)$ preserves the sparsity of their discretization. This paves the way for fast $O(N\ln N)$ and spectrally accurate mixed implicit-explicit time-marching of sets of linear and nonlinear equations in unbounded geometries.

Tables

Table 1. Table 1: Comparison of accuracy for the eigenmodes obtained by a spectral method on the infinite line ( S ​ B n 𝑆 subscript 𝐵 𝑛 SB_{n} or T ​ B n 𝑇 subscript 𝐵 𝑛 TB_{n} ), a Chebyshev collocation method, or first order finite differences with 256 modes or grid points. L o ​ p ​ t subscript 𝐿 𝑜 𝑝 𝑡 L_{opt} is the optimal mapping parameters for methods on the infinite line or the optimal size of the domain for the two other methods. N < 0.01 subscript 𝑁 absent 0.01 N_{<0.01} denotes the number of eigenvalues accurate within a relative error of 10 − 2 superscript 10 2 10^{-2} ; Q opt = N < 0.01 / 256 subscript 𝑄 opt subscript 𝑁 absent 0.01 256 Q_{\mathrm{opt}}=N_{<0.01}/256 is the corresponding proportion of accurate eigenvalues.

Method	$L_{opt}$		$Q_{\mathrm{opt}}$ (%)
$S{B}_n\left(y\right)$	$11$	$132$	$51.6\%$
$T{B}_n\left(y\right)$	$11$	$132$	$51.6\%$
Collocation $T_n\left(y\right)$	$11.75$	$75$	$29.3\%$
Finite Differences	$14$	$7$	$2.7\%$

Equations

$$\mathcal{M}{\partial }_t\mathbf{\psi }=\mathcal{L}\mathbf{\psi }+\mathcal{N}(\mathbf{\psi })$$

$${\mathcal{B}}_m\mathbf{\psi }(\mathbf{r})=b_m(\mathbf{r})\mathrm{for}\mathbf{r}\in \partial \Omega .$$

$${\mathbf{\psi }}_i(\mathbf{r})=\sum _{k=1}^N{\tilde{\mathsf{X}}}_{i,k}{\varphi }_k(\mathbf{r}).$$

$$\text{mathbfsf}M{\partial }_t\tilde{\mathsf{X}}=\text{mathbfsf}L\tilde{\mathsf{X}}+\mathsf{N}[\tilde{\mathsf{X}}].$$

$$\lambda \text{mathbfsf}M{\tilde{\mathsf{X}}}_a=\text{mathbfsf}L{\tilde{\mathsf{X}}}_a.$$

$$\left(\text{mathbfsf}M-\tau \text{mathbfsf}L\right){\tilde{\mathsf{X}}}^{(n+1)}=\text{mathbfsf}M{\tilde{\mathsf{X}}}^{(n)},$$

$$⟨f,g⟩=\frac{1}{\pi }{\int }_0^{\pi }f(\theta )g(\theta )\mathrm{d}\theta ,$$

$$f(\theta )=\sum _{n\ge 0}{\tilde{f}}_n\cos (n\theta ),$$

$$\mathrm{for}n=0:{\tilde{f}}_0$$

$$\mathrm{for}n>0:{\tilde{f}}_n$$

$$f(\theta )=\sum _{n\ge 1}{\hat{f}}_n\sin (n\theta ),$$

$${\hat{f}}_n$$

$$\frac{⟨\cos (m\theta ),\mathcal{L}\cos (n\theta )⟩}{\Vert \cos (m\theta )\Vert },\frac{⟨\sin (m\theta ),\mathcal{L}\cos (n\theta )⟩}{\Vert \sin (m\theta )\Vert },\frac{⟨\cos (m\theta ),\mathcal{L}\sin (n\theta )⟩}{\Vert \cos (m\theta )\Vert },\mathrm{or}\frac{⟨\sin (m\theta ),\mathcal{L}\sin (n\theta )⟩}{\Vert \sin (m\theta )\Vert },$$

$$\frac{⟨\sin (m\theta ),\mathcal{L}\sin (n\theta )⟩}{\Vert \sin (m\theta )\Vert }=\frac{1}{2}\left({\delta }_{n+1}^m+{\delta }_{n-1}^m\right)+\frac{m\left(1+{\left(-1\right)}^{m+n}\right)}{\pi (m-n+1)(m+n-1)}+\frac{m\left(1+{\left(-1\right)}^{m+n}\right)}{\pi (m-n-1)(m+n+1)}.$$

$$f(\theta )=f^{(s)}(\theta )+f^{(a)}(\theta ),$$

$$f^{(s)}(\theta )=f^{(s)}(\pi -\theta )\mathrm{and}{f}^{(a)}(\theta )=-f^{(a)}(\pi -\theta ).$$

$$f(\theta )=f^{(s)}(\theta )+f^{(a)}(\theta )$$

$$f(\theta )$$

$${\breve{\beta }}_n(\theta )=\{\begin{array}{rlr}& \cos (2p\theta )& \mathrm{for}n=2p+1,\\ & \sin (2p\theta )& \mathrm{for}n=2p,\end{array}$$

$${\mathring{\beta }}_n(\theta )=\{\begin{array}{rlr}& \cos ([2p-1]\theta )& \mathrm{for}n=2p-1,\\ & \sin ([2p-1]\theta )& \mathrm{for}n=2p,\end{array}$$

$$\displaystyle\Big{\{}\breve{\beta}_{n}(\theta)\Big{\}}_{n\geq 1}=\Big{\{}1,\sin 2\theta,\cos 2\theta,\sin 4\theta,\cos 4\theta,\cdots\Big{\}}\,,$$

$$\displaystyle\Big{\{}\mathring{\beta}_{n}(\theta)\Big{\}}_{n\geq 1}=\Big{\{}\cos\theta,\sin\theta,\cos 3\theta,\sin 3\theta,\cos 5\theta,\sin 5\theta,\cdots\Big{\}}\,.$$

$$\frac{⟨{\mathring{\beta }}_m(\theta ),\mathcal{L}{\breve{\beta }}_n(\theta )⟩}{\Vert {\mathring{\beta }}_m(\theta )\Vert },\frac{⟨{\breve{\beta }}_m(\theta ),\mathcal{L}{\mathring{\beta }}_n(\theta )⟩}{\Vert {\breve{\beta }}_m(\theta )\Vert }.$$

$$\mathcal{L}=\sum _{i=1}^{N_{\mathcal{L}}}\left({\alpha }_i\cos (n_i\theta )+{\beta }_i\sin (m_i\theta )\right){\left(\frac{\mathrm{d}}{\mathrm{d}\theta }\right)}^{{\gamma }_i}$$

$$y=L\cot \theta ,$$

$$\theta =\mathrm{arccot}\left(y/L\right),$$

$${\partial }_y=-\frac{{\sin }^2\theta }{L^2}{\partial }_{\theta }.$$

$$\psi (y)=\psi (y(\theta ))=f(\theta ).$$

$$T{B}_n(y)=\cos \left(n\mathrm{arccot}\left(y/L\right)\right),$$

$$S{B}_n(y)=\sin \left(n\mathrm{arccot}\left(y/L\right)\right).$$

$$\psi (y)\approx \sum _{n=0}^N{\tilde{\psi }}_{n}T{B}_n\left(y\right),$$

$$\psi (y)\approx \sum _{n=1}^N{\hat{\psi }}_{n}S{B}_n\left(y\right).$$

$${\int }_0^{\pi }f(\theta )\mathrm{d}\theta =\frac{1}{L}{\int }_{-\infty }^{\infty }\frac{\psi (y)}{1+y^2}\mathrm{d}y$$

Full text

Spectral methods; Sparse solvers; Chebyshev functions; Unbounded domain; Localized dynamics

I Introduction

Different numerical strategies (see boyd for a review) are available to solve differential equations on an infinite interval $y\in(-\infty,\infty)$. The first subcategory of methods consists in artificially introducing a bounded interval, the width of which is large compared to the extension of the exact solution. Over this finite interval, the classic weaponry of finite differences methods or spectral methods boyd ; fornberg ; orszag may be deployed. Other approaches rely on a representation of the solution using a basis of functions of infinite support. Usual candidates contain the Hermite functions, the Whittaker sinc function, or rational Chebyshev functions. In this paper, the focus is on the latter rational Chebyshev functions since fast Fourier transform algorithms can be advantageously utilized to compute the direct and inverse transforms from physical to spectral space.

Fast transform algorithms render amenable mixed implicit-explicit (IMEX) time-marching schemes for coupled systems of differential equations

[TABLE]

with linear operators $\mathscr{M},\mathscr{L}$ and where a popular choice is to treat pseudo-spectrally the nonlinear terms $\mathscr{N}$. Often, the system of differential equations are accompanied with a set of $M_{b}$ linear boundary conditions, labelled by the integer $m\in[1,M_{b}]$:

[TABLE]

The stability and accuracy of implicit time-stepping schemes render them very desirable over explicit schemes, but comes at the cost of one or several linear algebraic solves of the form $\mathbf{Ax}=\mathbf{b}$ for each time-step. As the resolution $N$ increases, such solves become prohibitively costly and necessitate $O(N^{3})$ operations if the chosen discretization results in a dense representation of $\mathscr{L}$. An implicit treatment of the linear terms, as in the popular IMEX schemes, hence advocates for a spectral decomposition that yields a sparse representation of both the operators $\mathscr{L}$ and $\mathscr{M}$. Additionally, solving the corresponding generalized eigenproblem $\lambda\mathscr{M}\boldsymbol{\psi}=\mathscr{L}\boldsymbol{\psi}$ and (2) in order to analyze the linear dynamics also benefit from sparse, and ideally well-conditioned, discretized operators.

In the case of a scalar differential equation, Boyd boydJCP87 has shown that the rational Chebyshev functions potentially provide a banded representation for operators $\mathscr{L}$ composed of multiplications by monomials and differential operators $y^{a}\left(\partial/\partial y\right)^{b}$. However, sparsity relies on the parity of $\mathscr{L}$, and is usually ruined if $\mathscr{L}$ is composed of operators of different parities.

In this paper, we introduce hybrid Chebyshev expansions based upon the definition of two families of function $TS_{n}(y)$ and $ST_{n}(y)$. These are obtained by interweaving Chebyshev rational functions $TB_{n}(y)$ and $SB_{n}$. Upon leveraging these hybrid expansions, we expose a method that yields sparse representations for rational differential operators, even in cases where parity mixing would destroy the sparsity with regular Chebyshev expansions. Following the spirit of this approach for a single scalar equation, a recipe for obtaining a sparse discretized system for sets of coupled differential equations is also exposed.

For the sake of completeness and clarity, the authors have elected to use a pedagogical tone. As such the first sections will be familiar to seasoned users of spectral methods. The complexity of the examples developed is gradually ramped up throughout the paper. The role of these examples is hence to illustrate the mechanics of our expansions in different situations, so that to provide the reader with models that he can adapt to his own problem.

The paper is organized as follows. A brief summary of spectral methods and discretization is presented in section II. The rational Chebyshev functions and hybrid expansions are presented in section III. Their convergence properties are established and their utility for producing sparse discretized system is illustrated on a fundamental textbook example: the quantum harmonic oscillator. Section IV is dedicated to the resolution of scalar equations with mixed-parity operator, through the example of the anharmonic oscillator. The method is extended to sets of coupled equations in Section V, through the geophysical fluid dynamics example of the equatorial $\beta-$plane equations. Section VI is dedicated to the generalisation to multidimensional domains. The applicability of the method to time-stepping schemes for non linear equations is described in section VII and illustrated on the example of the Kelvin-Helmholtz instability in shear flows. Finally, section VIII is composed of concluding remarks. Throughout, details of numerical approach are relegated to the appendices.

II From equations to discrete representations

Three classes of spatial discretization methods are generally available for governing equations (1) and (2). The venerable finite differences method, first introduced in 1910 by Richardson richardsonPRSA10 evaluates derivatives on a mesh and can be extended to the finite volume method on an irregular mesh leveque . Finite spectral element methods pateraJCP84 decompose $\Omega$ in sub-domains where the unknowns are approximated by a small sample of basis functions. Finally, spectral methods, of interest here, express the unknowns as truncated expansion of $N$ orthogonal basis functions \big{\{}\phi_{k}(\boldsymbol{r})\big{\}}_{1\leq k\leq N} on the whole domain $\Omega$ boyd ; fornberg ; orszag :

[TABLE]

Upon substitution of the truncated expansion above into the governing equations (1), projecting each individual equation onto the \big{\{}\phi_{j}\big{\}} basis yields the following algebraic system:

[TABLE]

In the context of spectral methods boundary conditions can be enforced following three routes. Galerkin methods use basis functions $\phi_{k}(\boldsymbol{r})$, possibly tailored by basis recombination, that intrinsically obey the boundary conditions on $\partial\Omega$, so that no additional step is required for the whole expansion to comply with boundary conditions. For basis functions that violate the boundary conditions, two approaches are available. The Collocation method explicitly solves the governing equations on the physical space, i.e. on the Gauss-Lobatto or Gauss-Chebyshev grids, and enforces the boundary condition on the boundary locations $\partial\Omega$. The tau method evaluates the boundary conditions in coefficient space and incorporates the resulting tau lines in the matrices $\mathbfsf{M}$ and $\mathbfsf{L}$ at the expense of higher order projections of the governing equations. A more detailed discussion is dedicated to boundary conditions in section VI.

The linear limit of the algebraic system (4) is obtained by setting $\mathsf{N}[\widetilde{\mathsf{X}}]=0$. In the case of autonomous sets of equations, a generalized eigenvalue problem where solutions are sought in the form of normal modes $\widetilde{\mathsf{X}}(t)=\widetilde{\mathsf{X}}_{a}\exp(\lambda t)$ with growth rate $\lambda$ can be solved to characterize the linear dynamics:

[TABLE]

When the full spectrum is sought after, QR/QZ methods are available in linear algebra libraries (e.g. LAPACK). Regardless of the sparsity of the matrices, such methods are very costly with an $O(N^{3})$ operations count. Moreover, it is often the case that the growthrate and frequency of the most unstable mode only is of primary interest. If this is so, using sparse iterative methods such as Arnoldi’s method will prove much cheaper and faster, typically of complexity $O(N)$ for sparse matrices. Iterative methods are therefore appealing for continuation in parameters. Granted, these methods suffer a certain lack of robustness: they usually require to be initialized with a good guess and commonly miss eigenvalues. However, a possible approach (previously described in boyd ) for a practical problem such as computing marginal stability curves or surfaces, or computing the optimal growthrate in a parameter space, is to first sweep through the parameter space along a coarse mesh and to compute the full spectra at a modest resolution using a QR/QZ algorithm. Once a first intuition is gained, a much finer sweeping can be carried out using iterative methods and a higher resolution (potentially ruling out under resolved modes, for instance).

An alternate approach suitable to the study of instabilities consists in time-stepping the linearized equation: one may then observe the growth of the most unstable mode. A large variety of time-stepping schemes are now available. Having applications in fluid mechanics in mind, the focus in this paper is on developing, for systems on an infinite line, spectral methods that are suitable for IMEX time-steppers. IMEX schemes evaluate the linear term in (4) implicitly at the end of the time-step, as opposed to evaluation at the start of the time-step for explicit time-stepping schemes. This result in a greater numerical stability and accuracy in comparison with their explicit counterparts, particularly desirable when stiff operators are present. As a pedagogical example, the simplest implicit time-stepper is the first order accurate Backward Euler scheme, which obtains the value of the unknown $\widetilde{\mathsf{X}}^{(n+1)}$ after a time-step $\tau$ by solving:

[TABLE]

where the value at the beginning of the time-step $\widetilde{\mathsf{X}}^{(n)}$ is known. Higher-order L-stable implicit schemes ascherANM97 ; groomsJCP11 , and more recently memory-efficient schemes cavaglieriJCP15 have been developed. The performance bottleneck of such schemes is the complexity of the linear solves of the form given in equation 6. Such linear algebraic systems are costly to solve when the matrix is dense: the complexity of the solve algorithm is of order $O(N^{3})$. In strong contrast, for a discretization that yields a banded operator $\left(\mathbfsf{M}-\tau\mathbfsf{L}\right)$, this complexity is drastically reduced to $O(N)$. From this observation, a successful implementation of an implicit algorithm is conditional on the existence of sparse matrix representations $\mathbfsf{M}$ and $\mathbfsf{L}$ for operators $\mathscr{M}$ and $\mathscr{L}$.

Beyond the linear dynamics, one may investigate the fully nonlinear dynamics, for instance in view of understanding the saturation of an instability. To avoid a costly direct computation of the non linear term in spectral space by means of a convolution product, the pseudo-spectral approach is utilized where potential derivatives are first evaluated in spectral space and products are then computed in physical space. The complexity of this procedure is dominated by the complexity of the transforms required between the physical and spectral space. When implemented as matrix multiplications, these transforms usually necessitate $O(N^{2})$ operations. This is largely the reason why spectral methods did not gain popularity until a fast Fourier transform (FFT) algorithm was published in 1965 by Cooley and Tukey cooleyMC65 , revisiting an idea of Gauss gauss . The FFT algorithm implements transforms in $O(N\ln N)$ operations only.

The great strength of spectral methods is their accuracy: the truncated expansion employed in spectral methods converges exponentially towards the true solution, when analytic, as opposed to algebraically with finite differences. Further, spectral methods have negligible spurious numerical dissipation, a well known feature of finite differences methods. A conservative numerical method is highly desirable for weakly dissipative equations, such as high Reynolds number fluid dynamics problems. The accuracy of spectral method comes at a price. First, they suffer from a poor flexibility concerning the geometry of the domain. Regular domains (e.g. rectangles, parallelepipeds, spheres, cylinders, etc.) are good candidates for a spectral treatment, but spectral methods are ill-suited for irregular geometries. Second, we already mentioned the necessity of transforms: this point strongly advocates for Chebyshev or Fourier basis, which both utilize fast Fourier transform algorithms. Third, the discretized representation of the operators are often dense.

As a conclusion of these observations, the emphasis of this paper is on choices of rational Chebyshev function expansions for functions on the infinite line. The central idea here is that sensibly chosen hybrid expansions obtained by interleaving different bases preserve the sparsity of discretized operators and exploit existing FFT algorithms.

III The rational Chebyshev functions

In this paper, we introduce the basis of rational Chebyshev functions as remapped cosine and sine functions on the infinite line. Not only does this approach lead to an insightful definition of the collocation points on the infinite line, but it also provides an understanding of which operators will lend themselves to a sparse discrete representation when approximated under the basis of rational Chebyshev functions. We thus begin by focussing on functions on the interval $[0,\pi]$ prior to the remapping.

III.1 Trigonometric basis functions on a finite interval

We consider $\mathcal{C}^{0}\left(\left[0,\pi\right]\right)$, the Hilbert space of real valued continuous functions on the interval $[0,\pi]$, that are also piecewise of class $\mathscr{C}^{1}$ (i.e. with a continuous first derivative) and that vanishes at the end of the interval $f(0)=f(\pi)=0$. We define the inner product:

[TABLE]

and the corresponding 2-norm $\left\|f\right\|=\sqrt{\left\langle f,\,f\right\rangle}$. The families of functions $\left\{\cos(n\theta)\right\}_{n\in\mathbb{N}}$, indexed by nonnegative integers, and $\left\{\sin(n\theta)\right\}_{n\in\mathbb{N}^{*}}$, indexed by positive integers $\mathbb{N}^{*}$, each form an orthogonal and countable basis of $\mathcal{C}^{0}\left(\left[0,\pi\right]\right)$. Hence a function $f\in\mathcal{C}^{0}\left(\left[0,\pi\right]\right)$ can be expressed as a cosine series:

[TABLE]

with

[TABLE]

or alternatively as a sine series:

[TABLE]

with

[TABLE]

For functions in $\mathcal{C}^{0}\left(\left[0,\pi\right]\right)$, both expansions converge uniformly to the function $f$. This does not imply that the expansions are equivalent: a discussion of the convergence properties of different representations, based on the analyticity of the possible extensions of $f(\theta)$ to the interval $[-\pi,\pi]$ (e.g., see  boyd ; fornberg ; orszag ), can be found in section III.5 below.

III.2 Hybrid trigonometric bases

Pure cosine or sine bases (presented above in equations (8,10)) suffer the curse of dense discretization for symmetry breaking operators. For instance, consider $\mathscr{L}=\cos\theta+\sin\theta$ as example of an operator with no particular symmetry around $\pi/2$. To avoid confusion with the symmetry around $\theta=0$ that stems in our discussion about convergence (section III.5 and appendix X), we emphasize that the symmetry that pertains to the present discussion is the one around $\pi/2$, the center of the segment $[0,\pi]$. Regardless of the selection of the pure basis (8) or (10) for expansion and projection, we would invariably obtain dense representations:

[TABLE]

all of which are band-unlimited. For example, the latter element of this list, represented on figure 1, yields:

[TABLE]

To circumvent the apparition of dense matrices, the strategy is to introduce mixed trigonometric bases vasilJCP08 , described below. The core idea to this paper then consists in generating the corresponding mixed Chebyshev rational function bases when working on an unbounded interval.

The continuous functions $f\in\mathcal{C}^{0}\left(\left[0,\pi\right]\right)$ may be uniquely decomposed into symmetric and antisymmetric components about the center of the segment $[0,\pi]$, namely $\theta=\pi/2$:

[TABLE]

In contrast to the “pure” expansion in sine and cosine presented above in equations (8,10), we define the two “hybrid” (or “mixed”) cosine-sine expansions:

[TABLE]

As deduced from these mixed representations, each component $f^{(s)}$ and $f^{(a)}$ has two possible basis representations. For instance, $\left\{\cos\left(2p\theta\right)\right\}_{p\in\mathbb{N}}$ and $\left\{\sin\left(\left[2p+1\right]\theta\right)\right\}_{p\in\mathbb{N}}$ each form an independent basis for symmetric functions about $\pi/2$, whereas, $\left\{\cos\left([2p+1]\theta\right)\right\}_{p\in\mathbb{N}}$ and $\left\{\sin\left(2p\theta\right)\right\}_{p\in\mathbb{N}^{*}}$ form an independent basis for antisymmetric functions about $\pi/2$. Since this paper leverages the virtue of mixed expansions (15a,15b) and makes a heavy use of them in producing sparse discretizations, each of these composite expansions is seen hereafter as a unique expansion by introducing the following compact notations:

[TABLE]

The two distinct hybrid bases obtained, $\{\breve{\beta}_{n}\}$ and $\{\mathring{\beta}_{n}\}$, can be thought of as interleaved collections of sines and cosines:

[TABLE]

or more explicitly:

[TABLE]

Both bases are orthogonal sets, as proven in appendix XI. The benefits of using such functions become evident by revisiting the introductory example of this paragraph, $\mathscr{L}=\cos\theta+\sin\theta$. From the trigonometric identities for products, we deduce easily that the following matrix representations for $\mathscr{L}$ (one of which is represented on figure 1) are both sparse:

[TABLE]

Therefore, a sparse discretization for $\mathscr{L}$ is obtained by considering an expansion in $\{\breve{\beta}_{n}\}$ and a projection on the $\{\mathring{\beta}_{n}\}$ basis (or vice versa). This observation can be generalized in the following proposition.

Proposition 1 –

Let $\mathscr{L}$ be a linear operator of the form

[TABLE]

where $N_{\mathscr{L}}\in\mathbb{N}$, and for all integers $1<i<N$: $\alpha_{i},\beta_{i}\in\mathbb{R}$, $\gamma_{i}\in\mathbb{N}$, and $(n_{i},m_{i})\in\mathbb{N}^{2}$ such that all integers $n_{i}$ and $m_{i}$ are of the same parity (even or odd). If all these integers are:

• •

odd, then $\frac{\left\langle\mathring{\beta}_{m}(\theta),\,\mathscr{L}\breve{\beta}_{n}(\theta)\right\rangle}{\left\|\mathring{\beta}_{m}(\theta)\right\|}$ and $\frac{\left\langle\breve{\beta}_{m}(\theta),\,\mathscr{L}\mathring{\beta}_{n}(\theta)\right\rangle}{\left\|\breve{\beta}_{m}(\theta)\right\|}$ are band-limited matrices,

• •

even, then $\frac{\left\langle\breve{\beta}_{m}(\theta),\,\mathscr{L}\breve{\beta}_{n}(\theta)\right\rangle}{\left\|\breve{\beta}_{m}(\theta)\right\|}$ and $\frac{\left\langle\mathring{\beta}_{m}(\theta),\,\mathscr{L}\mathring{\beta}_{n}(\theta)\right\rangle}{\left\|\mathring{\beta}_{m}(\theta)\right\|}$ are band-limited matrices.

Proof–

These matrices are easily built as linear combinations of matrix product of the banded elementary matrices $\mathcal{S}$, $\mathcal{C}$, and $\mathcal{D}$ given in appendix XII.1.

The requirement that the $n_{i}$ and $m_{i}$ be of given parity will be naturally satisfied by operators originating from a remapping of rational operators on the infinite line, as discussed in the following sections.

For completeness, we close up this section by mentioning alternative expansions: analogous to the $\{\breve{\beta}_{n}\}_{n\geq 1}$ basis is the ,$\{\exp\left(\mathrm{i}[2p+1]\theta\right)\}_{p\in\mathbb{Z}}$ basis. Similarly, the $\{\exp\left(\mathrm{i}[2p+1]\theta\right)\}_{p\in\mathbb{Z}}$ basis is analogous to the $\{\mathring{\beta}_{n}\}_{n\geq 1}$ basis. The changes of basis are straightforward, so that convergence properties (discussed below in section III.5) will be unaffected and choosing between the trigonometric and the complex exponential bases is purely a matter of taste. In the rest of this paper, we discuss and utilize the hybrid trigonometric expansions $\{\breve{\beta}\}$ and $\{\mathring{\beta}\}$ only.

III.3 Rational Chebyshev basis functions on the infinite line

Given the mixed representations on the $[0,\pi]$ interval, following Cain et al. cainJCP84 we now map this interval onto the infinite line using:

[TABLE]

where $L>0$ is the mapping parameter. Derivatives on the two intervals are connected according to

[TABLE]

By use of composite functions, a function on the infinite line $\psi(y)$ can be bijectively mapped on a function $f(\theta)$ on the segment $\theta\in[0,\pi]$:

[TABLE]

Hence, we define the Chebyshev functions $TB_{n}(y)$ and $SB_{n}(y)$ as the mapped cosine and sine functions stretched on the infinite line:

[TABLE]

From their definition and our discussion above, the Chebyshev functions inherit the properties of the trigonometric functions from which they are derived: they form two independent orthogonal bases on the infinite line and a truncated sum can approximate an arbitrary continuous bounded function:

[TABLE]

These expansions are the infinite line counterparts of the cosine and sine expansion given in equations (8,10). The $TB_{n}$ (respectively, $SB_{n}$) coefficients $\widetilde{\psi}_{n}$ (respectively, $\widehat{\psi}_{n}$) are equal to the cosine (respectively, sine) coefficients $\widetilde{f}_{n}$ (respectively, $\widehat{f}_{n}$). However, they can be computed directly on the infinite line, without mapping back to the interval $[0,\pi]$, by using

[TABLE]

to generalize equations (9,11) into:

[TABLE]

III.4 Hybrid Chebyshev bases

By means of the mapping (21), our partition of symmetric and asymmetric functions around $\theta=\pi/2$ on the interval $[0,\pi]$ results in a partition of functions $\psi(y)$ into odd and even components $\psi^{(e)}(y)$ and $\psi^{(o)}(y)$ around $y=0$. We emphasize here that in the following, symmetry properties are implicitly relative to $\theta=\pi/2$ for functions on the $[0,\pi]$ segment, and relative to $y=0$ on the infinite line $(-\infty,\infty)$. The odd and even components of $\psi(y)$ can be expanded independently with either a $TB_{n}$ or $SB_{n}$ basis. Hence, the analogous extension of the hybrid sine-cosine expansions (eq. 15) are the hybrid $SB$-$TB$ expansions:

[TABLE]

Following section III.2, we introduce the following notations for the two bases $ST(y)$ and $TS(y)$:

[TABLE]

More explicitly,

[TABLE]

such that the expansions (15a,15b) become respectively:

[TABLE]

In the remainder of the paper, we discuss how the sparsity of discretized systems is affected by the choice of this expansion. Particularly, we highlight cases when a mixed expansion is necessary to conserve sparsity, by extending our Proposition 1 that pertains to $\theta\in[0,\pi]$ onto the infinite line $y\in(-\infty,\infty)$ by means of the mapping $y=L\cot\theta$.

III.5 Choosing between $\{TB_{n}\}$, $\{SB_{n}\}$, $\{TS_{n}\}$, $\{ST_{n}\}$ expansions and alternatives

At this stage, it is necessary to discuss the global relevance of hybrid $\{TS_{n}\}$ and $\{ST_{n}\}$ expansions for a given function $\psi(y)$. Granted the fact that they prove useful in providing sparse discretization of parity-flipping operators, the question arises of under which conditions they inherit the convergence properties of the regular $\{TB_{n}\}$ and $\{SB_{n}\}$ expansions? The answer to this question is related to the analyticity of the continuations of the remapped function $f(\theta)=\psi(L\cot\theta)$ to the $[-\pi,\pi]$ interval. A detailled discussion is given in appendix X, and boils down to the following conclusions.

Exponentially decaying functions that are dominated by $\exp(-b|y|^{a})$ with $a,b>0$ have subgeometrically converging coefficients for each of the four expansions. That is to say, hybrid expansions $\{TS_{n}\}$ and $\{ST_{n}\}$ will converge as fast as pure $\{TB_{n}\}$ and $\{SB_{n}\}$ expansions. For this reason, hybrid expansions should be given preference over the pure expansions in cases where they yield a sparse discretization, as discussed above and exemplified below. Therefore the question of using hybrid Chebyshev expansions instead of alternative methods such as expansions based on Hermite functions or Cardinal functions is really equivalent to the choice between regular Chebyshev expansions versus Hermite or Cardinal functions. Hence, we refer the reader to the abundant existing literature (a possible starting point is boyd ). The bottom line is that in most cases these three options are equivalent and none really outperforms the others from the perspective of an optimal representation. In practice however, as stated in our introduction, it is our opinion that Chebyshev rational functions have an edge whenever a sparse representation is available. Indeed, Hermite functions prove to be inpractical in the context of non linear dynamics due to the absence of implementation of fast transforms in common linear algebra libraries–despite the existence of such algorithms orszag86 ; boydJCP92 . Using the Cardinal functions collocation method does not require any transform but results in dense discretizations, which yield a much slower time-stepping or use of iterative eigensolvers.

III.6 The relevance of mapping before discretizing

As mentioned above, one could very well introduce an expansion of the forms given in equations (25) or (28) and discretize directly a system of the form:

[TABLE]

However, more often than not, the path of least resistance to obtain sparse matrices consists in first mapping the equation to the $[0,\pi]$ segment before proceeding to the discretization. The reason behind this is that our intuition with the action of operators on the trigonometric functions is much more developed than our insight with the Chebyshev functions. As we exemplify in the next paragraph with the case of the quantum harmonic oscillator, tweaking the equations with a view to obtaining a sparse system is much easier with trigonometric functions, whereas dealing directly with Chebyshev functions remains obscure.

III.7 A fundamental example: the Quantum Harmonic Oscillator

We illustrate the mechanics of the Chebyshev decomposition on the classic example of the Quantum Harmonic Oscillator (QHO). The QHO is a fundamental example of a particle propagating in a potential $V(y)$. The wave function $\psi$ and the eigenenergies $E$ of the particle are the solutions to the Schrödinger equation:

[TABLE]

The potential $V(y)=y^{2}$ is well-known to possess a point spectrum of exact analytic solutions, indexed by the integer $n\geq 0$:

[TABLE]

where $\mathrm{H}_{n}$ is the Hermite polynomial of degree $n$. As described above, we first solve this equation by mapping it to the segment $[0,\pi]$ using  (21). Doing so yields the corresponding eigenproblem for $f(\theta)=\psi(y)$:

[TABLE]

In this form, the last term with a $\sin^{2}\theta$ denominator leads to a dense discretization regardless of which expansion is adopted for $f(\theta)$. Therefore, we regularize the equation by the multiplication $\sin^{2}\theta$ throughout. Doing so readily yields the generalized eigenproblem of the form

[TABLE]

In this simple example, the discretization is a straightforward step: anticipating no particular difficulty, we adopt as a first guess a pure cosine expansion of $N$ terms for $f(\theta)$. By projecting the equations onto the $N$ basis functions $\left\{\cos(m\theta)\right\}_{0\leq m\leq(N-1)}$, we obtain a $N\times N$ system of the form:

[TABLE]

where the matrix coefficients $\mathsf{M}_{mn}$ and $\mathsf{L}_{mn}$ are given by (for $1\leq m,n\leq N$):

[TABLE]

The exact expression of matrices $\mathbfsf{M}$ and $\mathbfsf{L}$ is omitted, but was obtained by using the procedure given in appendix XII.1. Their banded form can be observed on figure 2. This sparse structure could have been easily anticipated, as the discretized operators $\mathbfsf{M}$ and $\mathbfsf{L}$ are composed exclusively by sparse representations on the cosines basis: i.e., multiplication by $\cos\theta$, $\cos^{2}(\theta)$, $\sin^{2}(\theta)$, second order derivation $\partial_{\theta\theta}$, and $\sin\theta\partial_{\theta}$. Note that an alternative choice consisting in expressing $f(\theta)$ as a sine expansion and projecting on the $N$ functions $\left\{\sin(m\theta)\right\}_{1\leq m\leq N}$ would yield a similarly sparse linear system $E\,\mathbfsf{M}\boldsymbol{\widehat{f}}=\mathbfsf{L}\boldsymbol{\widehat{f}}$ where:

[TABLE]

In figure 3, we compare different methods used to solved the QHO equations, all of which contain a free length parameter $L$: the problem is solved on the infinite line using either a $TB_{n}(y)$ expansion (cosines), or a $SB_{n}(y)$ expansion. For comparison, we also solve our problem on a large segment $y\in[-L,L]$ by means of either a classic Chebyshev collocation method or finite differences. A Dirichlet boundary condition $\psi(L)=\psi(-L)=0$ is imposed for the latter methods. In all cases, the number of grid points (finite differences and collocation method) or modes (spectral methods) is set to $N=256$. The parameter $L$ is increased between $0.25$ and $20$ by increments of $0.25$. We display as a function of $L$ the number $N_{<0.01}$ of eigenvalues correct within a relative error of $0.01$ obtained by each method. For each method, an optimal value $L_{opt}$ is found which corresponds to a optimal ratio $Q_{opt}=N_{<0.01}/256$ of accurate eigenvalues. The results are summarized in Table 1. From a sheer accuracy criteria, the $SB_{n}$ and $TB_{n}$ functions on the infinite line are sensibly preferable to Chebyshev polynomials. All three spectral methods are substantially more accurate than the method of finite differences. Finally, note on figure 3 that the poor conditioning of differentiation matrices on the Gauss-Lobatto grid (the collocation grid for $T_{n}$ polynomials) is visible, as the condition numbers increases as $N^{4}$. Discretizations based on rational Chebyshev functions do not suffer this phenomenon and their condition number scales quadratically with $N$.

We emphasize that in the present example as well as in the following examples, the computation of the full spectrum using a standard QR/QZ algorithm and its comparison to the analytic spectrum had been intended to provide a purely pedagogical illustration of the accuracy of the method. For real life applications, the sparsity of discretizations using Chebyshev functions should be leveraged through the use of iterative eigensolvers or sparse solvers.

IV Parity splitting operators in a scalar equation

IV.1 Particle in a translated harmonic well: Chebyshev rational functions

We have illustrated how approximations using Chebyshev basis functions on the infinite line yields an efficient and accurate resolution of the QHO equation. However, the operators $y^{2}$ and $\partial_{yy}$ are both even and hence parity-preserving. In this section, we first illustrate how the sparsity of the approximation method is ruined in presence of both even and odd parity operators. A resolution is then provided to restore sparsity for a large class of polynomials and rational differential operators.

Our introductory example is the case of a particle in a asymmetric parabolic potential $V_{a}(y)=y^{2}-6y$. The potential $V_{a}$ is a mere translation of the harmonic potential accompanied by a renormalization of the ground state energy: applying the transformation $y=y^{\prime}+3$ and $E=E^{\prime}-9$, this problem can be mapped back to the QHO treated in the previous section. Hence we refer to this case as the “translated quantum harmonic oscillator” (TQHO) and we confront our approach against the known analytical solutions of this equation. Following a similar approach as above, we rewrite the Schrödinger equation for $\psi(y)$

[TABLE]

in $\theta$-space, and we regularize by multiplication by $\sin^{2}\theta$, yielding:

[TABLE]

We identify on this simple example the trouble that accompanies parity flipping operators. Equation (40) contains two classes of operators: the identity operator on the left-hand-side and the operator $\left(y^{2}-\partial_{yy}\right)$ preserve the parity of the function they act on. By contrast, the operator $-6y$ reverses the parity. As a consequence, were $\psi$ to be written as a “pure” $TB$ or $SB$ expansion, at least one term will be band-unlimited and the coupling matrix would be dense. Suppose we assume a pure cosine or sine expansion (8,10) for $f$, as we did for the QHO. A discrete set of equations of the form $E\,\mathbfsf{M}\boldsymbol{\widetilde{f}}=\mathbfsf{L}\boldsymbol{\widetilde{f}}$ could be obtained by projection onto the cosine or sine basis. However, by contrast with the QHO and regardless of which one of the four options we elect, this system would never yield a sparse representation. Indeed, the two possible representations of $\mathscr{L}$:

[TABLE]

are all dense matrices due the existence the parity-break term that results in a band-unlimited representation.

We note, however, that the differential operators present in equation 41 obey the requirement of our Proposition 1. Therefore, the sparsity can be restored if the odd and even components of $\psi$ are expanded by a different basis. A solution consists in assuming an hybrid expansion in sines and cosines of the form (15) for the symmetric and antisymmetric components $f^{(s)}$ and $f^{(a)}$. For instance, if we choose the basis $\{\breve{\beta}_{n}\}$ (i.e. a cosine and sine expansion (see equation 15a) for $f^{(s)}$ and $f^{(a)}$, respectively), then we can compute the matrices representing the action of the operators that compose $\mathscr{L}$ on this basis:

[TABLE]

all of which are sparse. With these building blocks in hand, the expression of which is given in appendix XII.1, the discretized system:

[TABLE]

is readily obtained and the sparse structure of matrices $\mathbfsf{L}$ and $\mathbfsf{M}$ is illustrated in figure 4:

[TABLE]

As can also be observed in figure 4 when compared to figure 2, neither the exponential convergence of the series nor the excellent accuracy of the method are adversely affected by our using of a mixed basis. This illustrates the discussion of convergence given in section III.5. For exponentially decaying solutions, expansions based upon the \Big{\{}\mathring{\beta}_{n}\Big{\}} basis or the \Big{\{}\breve{\beta}_{n}\Big{\}} basis yield a similarly sparse representations and converge equally rapidly.

IV.2 The generalized case

Having motivated, illustrated and benchmarked our method on a simple example, we now summarize and formalize the procedure to tackle a broader class of eigenproblems:

[TABLE]

where the operator $\mathscr{H}$ is a sum of rational fractions and derivatives of maximal order $N_{\mathscr{H}}$:

[TABLE]

and $P_{i}$ and $Q_{i}$ are polynomials. Note that even more general equations such as $E\mathscr{H}_{1}\psi=\mathscr{H}_{2}\psi$ (with $\mathscr{H}_{1}$ and $\mathscr{H}_{2}$ of the same form as given above in 47) can be considered. However, for the sake of brevity and clarity, we present the method in the case of equation (46) only.

The key idea is that eigenproblems of this form will necessarily fall back into the domain of validity of Proposition 1 after remapping (21). In contrast, we point out the following limitation to our technique: transcendental operators will not lend themselves to a sparse discretization. The method to obtain sparse representations for rational differential operators proceeds as follows.

As a first step, multiply the equation by $C(y)$, the least common multiple polynomial of the $Q_{i}(y)$ polynomials:

[TABLE]

where the polynomials $P^{\prime}_{i}(y)=C(y)P_{i}(y)/Q_{i}(y)$. As a second step, map the equation to $\theta$-space, using the mapping (21) and the chain rule (22). Determine the highest power $h_{p}$ of $\sin\theta$ present as a denominator: multiply the equation by $(\sin\theta)^{h_{p}}$ to cancel out all denominators, yielding to an equation of the form:

[TABLE]

Finally, obtain a sparse discretization of this system using a mixed expansion (15).

IV.3 Particle in an arbitrary asymmetric well

The general method described above is quickly illustrated on an example of greater complexity than the translated parabolic well considered in section IV.1. We consider a potential of the form:

[TABLE]

which corresponds to the Hamiltonian:

[TABLE]

of the form given in equation (47). We first eliminate the rational fractions:

[TABLE]

Mapping this equation to $\theta$-space and eliminating the $\sin\theta$ denominators by multiplying by $\sin^{4}\theta$ gives, in this case:

[TABLE]

This equation readily yields a sparse discretization when a mixed expansion is used for $f$. The sparse matrices involved are displayed in figure 5. The eigenfunctions of two different wells are computed to illustrate the variety of potentials amenable with our method.

V Coupled systems: the equatorial $\beta-$plane

V.1 The equatorial $\beta-$plane

We now motivate and illustrate the application of the method for coupled systems of equations, namely, the linearized equatorial $\beta-$plane model gillBOOK relevant for tropical atmospheric or oceanic dynamics. On the surface of a sphere covered with a shallow layer of fluid, the azimuthal velocity $u$, meridional velocity $v$ and height of the fluid layer $h$ obey the shallow water governing equations around the equator:

[TABLE]

where $x$ is the azimuthal coordinate, $y$ is the meridional coordinate, and $a$ is the Rayleigh friction coefficient.

As with the harmonic oscillator, we first present the analytical results against which the numerical method will be validated. In absence of friction, i.e., $a=0$, the solutions are waves and these equations can be solved analytically by considering normal modes $\left(U\left(x,y\right),V\left(x,y\right),H\left(x,y\right)\right)=\left(u(y),v(y),h(y)\right)\,\exp\left(\mathrm{i}\left[kx-\omega t\right]\right)$. The variables $u(y)$ and $h(y)$ can both be expressed as functions of $v(y)$:

[TABLE]

The variable $v$ itself is determined as the solution of a Weber equation, similarly to the quantum harmonic oscillator:

[TABLE]

As above, the solutions are parabolic cylinder functions and $\Omega$ is quantized:

[TABLE]

Each value $\Omega_{n}$ corresponds to three different modes of frequency $\omega_{n}^{(1)}$, $\omega_{n}^{(2)}$, and $\omega_{n}^{(3)}$ found as the roots of equation (56). The wave functions of each mode $u^{(1)}_{n}(y)$, $h^{(1)}_{n}(y)$, $u^{(2)}_{n}(y)$, $h^{(2)}_{n}(y)$, $u^{(3)}_{n}(y)$, and $h^{(3)}_{n}(y)$, are obtained by using formulas (55). For $n\neq 0,-1$ these modes represent (westward and eastward propagating) inertio-gravity waves and a (westward propagating) Rossby wave. Two exceptions appear to this picture: for $n=0$, one of the three roots, $\omega=-k$, must be discarded. Indeed, this root corresponds to an exponentially growing solution $v(y)=\exp\left(y^{2}/2\right)$. The remaining root represents the mixed Rossby-gravity (Yanai) wave. Further, the Kelvin waves are not captured by expression (57). Sometimes dubbed the “$n=-1$ mode”, these waves correspond to $v$ identically zero, yielding:

[TABLE]

V.2 Numerical treatment

When Rayleigh friction is present, the complex growth rate $\lambda$ of a normal mode of azimuthal wave number $k$ obeys the eigenproblem:

[TABLE]

Observe that variables form two groups: $u$ and $h$ form the first group and $v$ forms the second. The reason behind this segregation becomes clear when observing the blocks which form $\mathscr{L}$. In equation (59), the diagonal blocks (marked with rectangles) are composed of parity-preserving operators, whereas the off-diagonal blocks contain exclusively parity-flipping operators. Hence, in the spirit of the mixed expansions we have been using before in scalar equations, here we will use expansions of opposed parities in $y$ for different variables. As emphasized before, this remedy becomes clear after remapping the equation to $\theta$-space ($u^{\dagger}(\theta)=u(L\cot\theta)$, etc.):

[TABLE]

It becomes immediately clear that the sparsity of the discretization of the system is preserved if expansions of opposite parities are used for $(u^{\dagger},h^{\dagger})$ and $v^{\dagger}$. In the following, we present results obtained with a cosine expansion for $u^{\dagger}$ and $h^{\dagger}$ and a sine expansion for $v^{\dagger}$ (which yield to the structure depicted in figure 6 for matrices $\mathbfsf{M}$ and $\mathbfsf{L}$):

[TABLE]

We compare in figure 7 the frequencies obtained by solving the discretized eigenproblem to Rossby and inertio-gravity waves (57a), and Kelvin waves (58a). The accuracy of a spectral method based on Chebyshev functions yields an excellent agreement. We also display on figure 7 some eigenmodes. We illustrate how for a given $\Omega_{n}$, Rossby and inertio-gravity waves with different frequencies correspond to the roots of equation (56). These modes have a common azimuthal velocity $v$ but distinct meridional velocity $u$ and height of fluid $h$. We show plots of the amplitude of the profiles $|u(y)|$, $|v(y)|$, and $|h(y)|$ using a semi-logarithmic scale to illustrate the exponential decay of both variables as $|y|\rightarrow\infty$. Finally, the exponential convergence of the three spectra $|\widetilde{u}|$, $|\widetilde{v}|$, and $|\widetilde{h}|$ is shown on the last panel of figure 7.

We readily obtain a time-stepping scheme for the dynamics of equatorial waves by utilizing the matrices $\mathbfsf{M}$ and $\mathbfsf{L}$ with a implicit third-order Runge-Kutta scheme ascherANM97 . We show temporal dynamics in figure 8. A initial perturbation localized around a latitude $y_{0}$ travels across the equator and is reflected when reaching the opposite latitude $-y_{0}$. In absence of dissipation, the wavepacket then continues bouncing back and forth between the bounding latitudes $-y_{0}$ and $y_{0}$, illustrating a well known feature (see for instance gillBOOK ) of the equatorial region: this region behaves as a trapping wave guide.

VI Multidimensional domains

VI.1 The spirit of the method

The machinery of the Chebyshev functions on the infinite line is readily generalized for separable operators to bi- and tridimensional domains with at least one unbounded direction. For the sake of brevity, this paper describes only bidimensional domains $\left(x,y\right)\in\Omega$, the generalization to 3D being a mere bookkeeping exercise. As described in heinrichsMoC89 , the philosophy consists in decomposing variables as a truncated expansion of a well-suited family of basis functions $\alpha_{n}(x)$ and $\beta_{m}(y)$ along each direction in the following fashion:

[TABLE]

The procedure is as follows. First, select the bases $\alpha_{n}$ and $\beta_{m}$: Fourier modes should be used along periodic directions; Chebyshev polynomials are almost always the best choice along bounded directions; unbounded directions can be remapped to $[0,\pi]$, where cosines, sines, or the hybrid bases $\breve{\beta}$ or $\check{\beta}$ can be used, depending on the parity properties of the operators involved as discussed above.

Second, regularize the equations so that they contain only operators with a sparse representation. We have indicated a procedure for unbounded directions in section IV.2; for bounded directions where Chebyshev polynomials are used, follow the quasi-inverse method of  julienJCP06 equivalent to repeated integrations to remove all derivatives along these directions.

Third, decompose the separable operators into a sum of the form $\mathscr{L}=\sum_{i}\mathscr{P}_{i}\mathscr{Q}_{i}$, where the $\mathscr{P}_{i}$ (resp. $\mathscr{Q}_{i}$) operators act on the variable $x$ (resp. $y$) only. Compute the matrices that represent the action of these operators on the selected bases $\alpha_{n}(x)$ and $\beta_{m}(y)$:

[TABLE]

These sparse matrices form the building blocks for the discretization of $\mathscr{L}$ which may now be assembled using the Kronecker product to obtain the matrix representing $\mathscr{L}$:

[TABLE]

A similar treatment should be applied to the operator $\mathscr{M}$. The Kronecker product offers a simple way to compute the representation of operators $\mathscr{M}$ and $\mathscr{L}$ on the basis $\left(\alpha_{0}\beta_{0},\alpha_{0}\beta_{1},\dots,\alpha_{0}\beta_{N_{y}},\alpha_{1}\beta_{0},\dots\alpha_{N_{x}}\beta_{N_{y}}\right)$. The matrix of spectral coefficients $\widetilde{f}_{nm}$ can be vectorized into a long column vector of state

[TABLE]

which represents a stack of the columns of $\widetilde{f}_{nm}$. Finally, we obtain the familiar form $E\,\mathbfsf{M}\,\boldsymbol{f}=\mathbfsf{M}\,\boldsymbol{f}$.

We illustrate this method on the three cases that occur with at least one unbounded direction. In this section, we consider a 2D quantum harmonic oscillator on the infinite plane $\Omega=\mathbb{R}^{2}$ and on the infinite strip with finite width $\Omega=\left[x_{1},x_{2}\right]\times\mathbb{R}$. In the next section, the case of a domain with an infinite direction and a periodic direction (infinite cylinder) is considered through the example of the Kelvin-Helmholtz instability. This is investigated through a time-stepping method in contrast with the eigenvalue problems considered in the present section.

VI.2 Bidimensional benchmark: the 2D QHO

We consider a particle in the two-dimensional potential

[TABLE]

and obeying the 2D Schrödinger equation

[TABLE]

The motivation behind the choice of potential $V(x,y)$ is that this potential is nothing more than the simple potential $\bar{V}(X,Y)=X^{2}+Y^{2}/3$ rotated by an angle $\pi/3$ and translated by the quantity $(1,2)$. The resulting potential thus possesses no parity in $x$ nor $y$ and features a coupling term $xy$. Hence, equation (67) is a test problem of reasonable complexity for our method. Despite this apparent complexity, the results are easily checked against analytic solutions. Indeed, since translations and rotations are unitary transforms, the eigenvalues of $V$ are equal to the eigenvalues of $\bar{V}$, namely:

[TABLE]

Further, the eigenfunctions of $V$ are obtained by applying the same rotation and translation to those of $\bar{V}$:

[TABLE]

VI.3 The infinite plane $\Omega=\mathbb{R}^{2}$

On the infinite plane, since we have demonstrated the efficiency and accuracy of the Chebyshev approximation, we simply decompose our variable on two sets of Chebyshev functions along $x$ and $y$. As in the 1D case, more intuition is generally available when dealing with trigonometric function instead of the full-fledged Chebyshev functions: hence we remap both the coordinates $x$ and $y$ using:

[TABLE]

where $L_{x}$ and $L_{y}$ are two real mapping parameters that can be optimised for accuracy, depending on the problem. Our Schrödinger equation (67) for $\psi(x,y)$ is mapped into the following equation for $f(\chi,\theta)=\psi(x(\chi),y(\theta))$:

[TABLE]

after multiplication (regularization) by $\sin^{2}\chi\sin^{2}\theta$ to preserve sparsity. Due to the presence of parity mixing in both directions $\theta$ and $\chi$, we discretize this equation using mixed expansion in $\chi$ and $\theta$: for instance $\breve{\beta}(\chi)$ and $\breve{\beta}(\theta)$:

[TABLE]

Note that nothing would prevent us from using different expansions in each direction, for instance, $\breve{\beta}(\chi)$ and $\mathring{\beta}(\theta)$. Doing so, however, would unnecessarily demand to compute the action of the operators on two different bases, yielding twice as much matrix building compared to using a common basis for both $\chi$ and $\theta$. The discretized version of system (71), with $\boldsymbol{f}=\mathrm{vec}(\widetilde{f}_{nm})$:

[TABLE]

where matrices $\mathbfsf{A}$, $\mathbfsf{B}$, $\mathbfsf{C}$, $\mathbfsf{D}$, and $\mathbfsf{E}$ are defined in equations (43a,b) and (43c,d,e) and computed in appendix XII.1.

VI.4 The infinite strip

We now turn to solving equation (67) on an infinite strip $\Omega=[x_{1},x_{2}]\times\mathbb{R}$ with Dirichlet boundary conditions $\psi(x_{1})=\psi(x_{2})=0$. As usual, the infinite $y-$direction is remapped to $\theta\in[0,\pi]$ and the bounded $x-$direction is remapped to $\chi\in[-1,1]$ via:

[TABLE]

Since the spirit of this paper is to provide recipes for memory efficient discretization methods along the infinite line, it is desirable to couple them with memory efficient methods along the other directions too. When a truncated Chebyshev polynomial expansion is considered along a bounded direction, a sparse representation is obtained for operators containing only polynomials and derivatives. The quasi-inverse method orszag is employed that consists of integrating the governing equation or set of equations as many times as the highest derivative present. In our example, we will integrate twice to cancel out the second order $x$ derivative, so that equation (67) becomes:

[TABLE]

Observe that by integrating twice we have introduced two arbitrary integration constants $a$ and $b$. The value of these yet floating integration constants will be determined by imposing boundary conditions as the last step of our discretization. For the time being, we define the discrete representation of three useful operators in Chebyshev space, denoted $\mathbfsf{Q}$:

[TABLE]

These matrices are banded with respectively two-, three- and four-bands. We give their full expressions in the appendix XII.2 and show that given the discretized matrix $\mathbfsf{x}$ for $x$ we arrive at $\mathbfsf{Q}^{(1)}=\mathbfsf{Q}^{(0)}\mathbfsf{x}$ and $\mathbfsf{Q}^{(2)}=\mathbfsf{Q}^{(0)}\mathbfsf{x}^{2}$ with simply matrix multiplications. Similarly to the previous section, we easily obtain our matrices $\mathbfsf{M}$ and $\mathbfsf{L}$ by combining the sparse representations of operators along different directions with the Kronecker product:

[TABLE]

The last step consists in implementing boundary conditions: as mentioned above the first two coefficients of the expansion in $x$ are polluted by integration factors $a(\theta)$ and $b(\theta)$. As a consequence, the $2N_{y}$ lines in $\mathbfsf{M}$ and $\mathbfsf{L}$ that correspond to coefficients $f_{0m}$ and $f_{1m}$ (with $m\in[1,N_{y}]$) have to enforce Dirichlet boundary conditions, namely:

[TABLE]

which we write

[TABLE]

for a greater generality. Any boundary conditions can be implemented easily by using only matrix products in the following manner. We define a boundary-condition matrix $\mathbfsf{b}$ which is zero everywhere, except for the first two lines $b_{0n}$ and $b_{1n}$, and a diagonal matrix $\mathbfsf{1}^{(2)}$ equal to the identity everywhere except on the first two lines, which are identically zero. Define:

[TABLE]

The mechanics of this simple algebra, illustrated on figure 9, is the following: multiplying by $\mathbfsf{1}^{(2)}\otimes\mathbfsf{1}_{y}$ zeroes out the $2N_{y}$ lines of matrices $\mathbfsf{M}$ and $\mathbfsf{L}$ that correspond to projection onto the two lowest order Chebyshev modes that are polluted by arbitrary integration constants, as seen in equation (75). In $\mathbfsf{L}_{\mathrm{B.C.}}$, these lines are replaced with the boundary condition by adding $\mathbfsf{b}\otimes\mathbfsf{1}_{y}$. Eigenmodes for both the infinite strip and the infinite plane are shown on figure 10. The accuracy of this bidimensional example is omitted but has been checked to be as good as in the 1D case.

VII Nonlinear time-stepping: the Kelvin-Helmholtz instability

VII.1 Background parallel flow and formulation

We consider a streamwise invariant parallel flow $U_{0}(y)\boldsymbol{e}_{x}$ imposed by external forces in a periodic domain along the streamwise direction $x\in\left[0,2\pi\right[$ with an unbounded spanwise direction $y\in\mathbb{R}$. Perturbations $\boldsymbol{u}(x,y,t)=-\partial_{z}\psi\,\boldsymbol{e}_{x}+\partial_{x}\psi\,\boldsymbol{e}_{z}$ to the background flow obey the incompressible governing equation:

[TABLE]

Such flows are known to be potentially linearly unstable if they present an inflection point $\partial_{zz}U_{0}\neq 0$ (see for example drazin2004 ). The Reynolds number, measuring the strength of the viscous dissipation force, controls the stability of the shear flow, i.e., the flow stabilizes as $Re$ is lowered. In this section, the stability of the flow

[TABLE]

is analyzed by means of a time-stepping integration of the full nonlinear governing equation (82).

VII.2 Discretization

Following the usual prescription, we remap the equations to $\theta$-space and multiply out with the factor $\left(L^{4}\cos^{4}\theta+\sin^{4}\theta\right)^{3}$ to obtain:

[TABLE]

where $\mathscr{M}$ and $\mathscr{L}$ are linear operators and $\mathscr{N}(f)$ is the nonlinear advection term:

[TABLE]

We expand our variable $f$ as a Fourier series along the periodic direction $x$. Observing that parity mixing is present amongst the operators, we chose to represent $f$ as a hybrid cosine and sine series in $\theta$:

[TABLE]

As the operators $\mathscr{M}$ and $\mathscr{L}$ are separable, we easily obtain their discrete representation $\mathbfsf{M}$ and $\mathbfsf{L}$ by using Kronecker products. We compute the nonlinear term pseudospectrally by first evaluating derivatives $\partial_{x}f$, $\partial_{\theta}\nabla^{2}f$, etc., in spectral space. We then transform these variables to physical space where the products are computed, and transform back to spectral space. A $2/3$ dealiasing rule, appropriate for quadratic nonlinearities, is used. Additional details on the computation of $\mathbfsf{N}(f)$, particularly concerning the implementation of the transform between spectral and physical space for the functions $\breve{\beta}_{n}$ by means of FFT algorithms, are relegated to appendix LABEL:app:transforms.

VII.3 Semi-implicit time-stepping schemes

IMEX schemes are a popular class of time-stepping schemes in computational fluid mechanics: these schemes combine an implicit treatment of the linear terms, whereas nonlinearities are treated explicitly. This popularity stems from the salient features of the Navier-Stokes equations which typically contains stiff linear terms, advocating for the use of stable and accurate implicit schemes. An implicit treatment of the nonlinear term would be costly and ultimately requires iterative methods compared to direct solves. Fortunately an explicit treatment yields satisfactory results. The simplest IMEX scheme, used in the following, is the order one Backward Euler which obtains $f^{(i+1)}$ from $f^{(i)}$ after a timestep $\tau$ by solving the linear system:

[TABLE]

VII.4 Stability of the flow: numerical results

Equation (82) is investigated by time-integration starting from an initial condition consisting of small amplitude noise uniformly distributed over the resolved spatial scales. The time series of the kinetic energy $E_{K}$ is displayed on figure 11 for seven values of the Reynolds number $Re$. Three phases are distinguished.

During the initial phase, as most of the modes that compose the initial condition are stable, the total energy first plummets (for $t\lesssim 20$). During the second phase, an optimal mode emerges from the sea of decaying modes: this mode is either the most unstable (fastest growing) mode or the least damped mode, depending on the value of Re. The planform of the mode that dominates the temporal dynamics is depicted in figure 12 for $\mathrm{Re}=20$ and compared with the optimal mode obtained by means of the corresponding eigenproblem. As expected, the agreement is excellent. As long as the dynamics remains in the linear regime ($|\psi|\ll 1$). The total kinetic energy is dominated by the energy of the optimal mode alone which follows an exponential growth trend in time $E_{K}\propto\mathrm{exp}(2\lambda t)$. The optimal growth rate $\lambda(Re)$ can be fitted for each values of Re. Figure 11 shows how the optimal growth rate $\lambda(Re)$ varies as a function of $Re$. We obtain the critical Reynolds number $\mathrm{Re}_{c}\approx 11.066$ corresponding to the threshold of the instability by fitting the growthrate $\lambda$ with a polynomial in $\mathrm{Re}$ around the onset.

Finally, the instability saturates as $\psi$ becomes of order $O(1)$ (figure 11). As can be observed from the timeseries at modest supercriticality $\mathrm{Re}\leq 20$, the saturated state becomes steady after a short overshoot. The steady planform of the saturated state is shown in figure 12 for $\mathrm{Re}=20$ and displays a noticeable but vertically contained distortion compared to the linear optimal mode. By contrast, when the supercriticality is increase (e.g. for $\mathrm{Re}=40$), the saturated state becomes unsteady, showing intermittent bursts of activity, as reported in the bottom row of figure 12. The unstable layer of fluid around $y=0$ then sheds vortices to the quiescent surrounding regions. Here, we notice the limitations of our method. As no proper physical trapping mechanism is present, viscous damping remains the only source of control, thus the outwardly propagating vortices can violate the CFL criterion as they propagate fast enough to reach regions where the collocation grid becomes coarser and coarser. If viscosity acts on timescales longer than the time needed for these vortices to reach large values of $|y|$, these vortices will become underresolved before they are damped out. The accuracy and stability of the methods is then jeopardized. If no additional physics is introduced to resolve this issue, the solution is to increase the number of basis function along $y$ and to use a larger $L$ to increase the size of the region where the dynamics is safely resolved. Note that a similar increase in resolution would be necessary if one were to chose to introduce boundary conditions to tackle the problem of an unstable shear flow in an unbounded $y-$direction: these boundary conditions would have to be placed “far enough”, so that they do not affect the dynamics. Increasing the distance between boundaries would equally lead to increasing the number of grid points or Chebyshev polynomials in the $y-$direction.

VIII Conclusion

In this paper, we have illustrated on fundamental examples originated from topics as varied as quantum mechanics, geophysical fluid dynamics, or fluid mechanics, the implementation of sparse spectral methods for solving PDE systems in unbounded domains of arbitrary dimensions.

Central to our technique are hybrid expansions $\{TS_{n}\}$ and $\{ST_{n}\}$ (equation 28), composed of Chebyshev functions $\{TB_{n}\}$ and $\{SB_{n}\}$ which preserve the sparsity of the discretization of differential operators on the infinite line even in presence of parity flipping differential operators in the general form of polynomial or rational fractions. Chebyshev functions correspond to remapped $\cos$ and $\sin$ functions to the infinite line, so that our hybrid expansions are remapped bases of interleaved sines and cosines, $\{\breve{\beta}_{n}\}$ and $\{\mathring{\beta}_{n}\}$ defined in equations 17. For exponentially decaying functions, we have shown that our hybrid Chebyshev expansions inherit the exponential convergence properties of standard Chebyshev expansions.

Numerical analyses of systems of differential equations are greatly facilitated, or made possible at all, by the sparsity of the discretization of rational differential operators that results from this choice of hybrid bases. We have shown examples of investigation of the linear dynamics by solving sparse eigenproblems or by means of strongly accurate and stable implicit time-marching schemes, along with examples of analyses of non linear dynamics, tackled by means of popular IMEX schemes. The method is easily generalized to domains of higher dimensionality by means of Kronecker products.

The work in this paper was supported in part by the National Science Foundation under Grant DMS-1317666. The authors would like to acknowledge useful comments from Dr. Ian Grooms and two anonymous referees.

IX Bibliography

X Appendix: Convergence of the expansions.

We discuss in this appendix the speed of convergence of the expansion of functions $\psi(y)$ on the $TS(y)$, and $ST(y)$ bases (eq. 30). We first present tools to assess the speed of convergence of functions on the $[0,\pi]$ interval using the $\mathring{\beta}$ and $\breve{\beta}$ bases (eq. 18). We then discuss a few cases of functions on the infinite line, distinguished by their asymptotic behaviour at $|y|\rightarrow\infty$.

X.1 Convergence of Fourier series for $2\pi$-periodic functions.

It is well-known (see e.g. [21]) that a $2\pi$ periodic function which is of $\mathscr{C}^{k-1}$ class (i.e. with a continuous $k-1$ derivative) on $\mathbb{R}$ and piecewise of $\mathscr{C}^{k}$ class on $\mathbb{R}$ possesses Fourier coefficients $\tilde{f}_{n}=o(1/n^{k})$. This behaviour is referred to as algebraic convergence. By extension, functions of $\mathscr{C}^{\infty}$ class on $\mathbb{R}$ possesses Fourier coefficients that decay asymptotically faster than any power of $1/n$. Among the class of $\mathscr{C}^{\infty}$ functions on $\mathbb{R}$, the most coveted are analytic functions, the Fourier coefficients of which decay exponentially with $\tilde{f}_{n}\sim\exp(-an)$, with $a>0$. Such a convergence is dubbed geometric convergence. Remapping infinite intervals onto the $[0,\pi]$ segment shall typically bring us to consider an intermediate class of functions: non analytic $\mathscr{C}^{\infty}$ functions with essential singularities. For instance, $\exp(-\cot(\theta)^{2})$ is a $2\pi$-periodic function of class $\mathscr{C}^{\infty}$ on $\mathbb{R}$ with essential singularities at $\theta=n\pi$ (with $n\in\mathbb{Z}$). It has been shown in [22] that Fourier coefficients of $\mathscr{C}^{\infty}$ functions which possess a singularity of the form $\exp(-1/\left|\theta\right|^{k})$ decay like $\tilde{f}_{n}\sim\exp(-an^{2k/[2k+1]})$, with $a>0$.

X.2 Convergence of cosine, sine, or mixed trigonometric series on the $[0,\pi]$ segment.

Let $f(\theta)$ be a function of $\mathscr{C}^{\infty}$ class on $[0,\pi]$, that satisfies $f(0)=f(\pi)=0$. The behaviour of the trigonometric series of $f$ depends on symmetry considerations. To arm ourselves for manipulating such functions, we define the projectors on the symmetric and antisymmetric components of $f$ with respect to the center of the interval $\pi/2$ (see figure 13, panels a and c):

[TABLE]

Using these projectors, four distinct continuations of $f$ to the interval $[0,2\pi]$, with different symmetries, can be built. The four continuations (illustrated on figure 13, panels b and d) naturally coincide with $f$ on the $[0,\pi]$ segment and are as defined as follows for $\theta\in[-\pi,0]$:

[TABLE]

Each of these continuations corresponds to one of our four choices of basis: $\left\{\cos(2p\theta)\right\}$ or $\left\{\sin([2p+1]\theta)\right\}$ for the symmetric component $\mathscr{S}_{\pi/2}[f(\theta)]$, and $\left\{\cos([2p+1]\theta)\right\}$ or $\left\{\sin(2p\theta)\right\}$ for the antisymmetric component $\mathscr{A}_{\pi/2}[f(\theta)]$. Inspecting the parity around $\theta=0$, one easily identifies that $g_{1}$ corresponds to the pure cosine expansion, $g_{4}$ to the pure sine expansion, $g_{2}$ to the mixed $\breve{\beta}$ expansion, and $g_{3}$ to the $\mathring{\beta}$ expansion:

[TABLE]

The behaviour of these four expansions, and notably their speed of convergence, is bound with the class of the continued fonction $g_{i}$ they correspond to. We discuss several cases below.

Functions with essential singularities at [math] and $\pi$–

Functions with essential singularities at the end of the $[0,\pi]$ segment (common in the context of the remapped infinite line) are advantageous and give complete flexibility concerning the choice of the basis. Indeed, consider $f$ such that:

[TABLE]

Each of the for continuations $g_{1},\dots,g_{4}$ will be of class $\mathscr{C}^{\infty}$ on $\mathbb{R}$, with essential singularities. Therefore, the coefficients of each expansion (pure or mixed trigonometric) will converge at a subgeometric rate. In this case, the user has a total freedom of choice: mixed expansions will inherit the desirable properties of pure expansions.

Non singular functions–

Non singular functions are ironically the least flexible as far as the choice of the basis is concerned. Let the Taylor expansions of $f$ at $\theta=0$ and $\theta=\pi$ be:

[TABLE]

Our decomposition based on the symmetry around $\pi/2$ can be taken one step further, to distinguish between even and odd functions around [math]. The function $f$ can be decomposed into:

[TABLE]

with

[TABLE]

As one attempts to continue $f$ over $[-\pi,0]$ by building $g$ of the form:

[TABLE]

it becomes obvious that none of the four choices for $\sigma_{\mathscr{S}},\sigma_{\mathscr{A}}=\pm 1$ will yield an analytic function unless $\mathscr{S}[f(\theta)]$ and $\mathscr{A}[f(\theta)]$ are purely odd or purely even at $\theta=0$: some derivatives will be discontinuous at [math]. Such functions $f(\theta)$ will invariably have algebraically decaying Fourier coefficients, regardless of the expansion chosen.

Functions such that $\mathscr{S}[f(\theta)]$ and $\mathscr{A}[f(\theta)]$ have a definite parity around zero have potential for being continued into a $2\pi$-periodic analytic function on $\mathbb{R}$, but will offer little flexibility concerning the choice of basis that yields geometric convergence. As an illustration, consider $f$ such that $f(\theta)=\mathscr{SE}(\theta)+\mathscr{AO}(\theta)$. Building $g_{2}(\theta)=+\mathscr{S}_{\pi/2}[f(-\theta)]-\mathscr{A}_{\pi/2}[f(-\theta)]$ is the sole way to continue $f$ into an analytic function. As a consequence, a $\{\breve{\beta}\}$ expansion alone yields geometrically decaying coefficients.

X.3 Convergence of $SB$, $TB$, or mixed $TS$ or $ST$ bases on $(-\infty,\infty)$.

By virtue of the discussion above, the speed of convergence of pure $\{TB\}$ and $\{SB\}$ or hybrid $\{TS\}$ and $\{ST\}$ Chebyshev expansions for a function $\psi(y)$ is controlled by the asymptotic behaviour of $\psi$ as $|y|\rightarrow\infty$. Functions $\psi(y)$ dominated at $|y|\rightarrow\infty$ by an exponentially decaying function $\exp(-b\left|y\right|^{a})$ with $a,b>0$, are mapped onto functions $f(\theta)$ with essential singularities at both ends of the $[0,\pi]$ segment. Indeed, upon remapping $y=L\cot\theta$ and choosing $L=1$ for simplicity and without loss of generality, we find that $f(\theta)=\psi(y(\theta))$ is dominated around $\theta\downarrow 0$ by $\exp(-b/\left|\theta\right|^{a})$ (by virtue of $\cot\theta=1/\theta+O(\theta)$). In the light of our reasoning on the $[0,\pi]$ segment, the four mixed and hybrid expansions will converge subgeometrically.

Functions that decay algebraically $\psi(y)\sim 1/y^{k}$ for $|y|\rightarrow\infty$ will map to analytic functions on $[0,\pi]$. As stressed above, such a situation is unfortunate as it might not allow to choose freely the basis of expansion (unless one is to pay the price of algebraic convergence).

XI Appendix: Orthogonality of the hybrid bases $\{\breve{\beta}_{n}\}$ and $\{\mathring{\beta}_{n}\}$.

Recalling from equation (17a) that:

[TABLE]

the orthogonality of the basis is immediate. We prove that $\left\langle\breve{\beta}_{m},\,\breve{\beta}_{n}\right\rangle\propto\delta^{n}_{m}$, where $\delta^{n}_{m}$ is the kronecker delta symbol, by remarking the following: we have $\left\langle\breve{\beta}_{1},\,\breve{\beta}_{1}\right\rangle=1$, and for $p$ or $q$ not both equal to zero:

[TABLE]

Similarly, using equation (17b):

[TABLE]

one immediately gets the orthogonality $\left\langle\mathring{\beta}_{m},\,\mathring{\beta}_{n}\right\rangle=\frac{1}{2}\delta^{n}_{m}$:

[TABLE]

XII Appendix: Matrix elements

XII.1 Discretization along an unbounded direction: Chebyshev functions

This appendix documents a method to obtain discretized matrices for operators composed of superpositions of powers of cosines, sines, and derivatives. Such operators are typically obtained after the procedure described in IV.2 and exemplified in IV.3. We will consider the four matrices $\mathcal{S}^{sc}$, $\mathcal{S}^{cs}$, $\mathcal{C}^{cc}$ and $\mathcal{C}^{ss}$ as our elementary building blocks. These matrices represent the action of $\sin\theta$ on cosines and sines, and the action of $\cos\theta$ on cosines and sines, respectively. Their matrix elements are given by:

[TABLE]

From these building blocks, more complex operators are readily constructed using matrix products, with the only reservation that the lower right corner of the matrix, which corresponds to modes neighbouring the truncation, deserves specific attention. Indeed, should no care be taken, the lower right corner of the matrix product would be inexact, due to the proximity of truncation. A rule of thumb is that, when two matrices $\mathbfsf{L}_{1}$ and $\mathbfsf{L}_{2}$ represent the action of operators $\mathscr{L}_{1}$ and $\mathscr{L}_{2}$ and have a bandwidth $(1+2p)$ (with $p$ subdiagonals and $p$ superdiagonals), the action of the operator $\mathscr{L}_{3}=\mathscr{L}_{1}\mathscr{L}_{2}$ is not represented by the product $\mathbfsf{L}_{1}\mathbfsf{L}_{2}$, for the last $p$ columns of this product are spurious. This problem can be understood as some kind of aliasing related to linear operators and is easily circumvented by using a larger truncation $N^{\prime}>N$ to discretize $\mathscr{L}_{1}$ and $\mathscr{L}_{2}$, so that to obtain large square matrices $\mathbfsf{L}_{1}^{\prime}$ and $\mathbfsf{L}_{2}^{\prime}$ of size $N^{\prime}$. Then, one computes the product $\mathbfsf{L}_{3}^{\prime}=\mathbfsf{L}_{1}^{\prime}\mathbfsf{L}_{2}^{\prime}$ of size $N^{\prime}$, and finally one extracts a smaller square matrix $\mathbfsf{L}_{3}$ of size $N$ that rigorously discretize $\mathscr{L}_{3}$. Below we denote with an overbar $\overline{\mathbfsf{L}}^{N}$ the extraction from a given square matrix $\mathbfsf{L}$ of a square submatrix of size $N$. Using our prescription, we can safely use a matrix product to obtain the matrices $\mathbfsf{A}$, $\mathbfsf{B}$, $\mathbfsf{C}$, $\mathbfsf{D}$ and $\mathbfsf{E}$, defined in equation (43), as products of elementary matrices $\mathcal{E}$ of suitably enlarged size $N^{\prime}>N$. We distinguish below the case of parity-conserving and parity-flipping operators.

The case of parity-conserving operators

Our prescription is first illustrated with the parity conserving operator $\mathscr{D}=\sin^{2}\theta$, in view of obtaining its discretization $\mathbfsf{D}$ on the basis $\{\breve{\beta}_{n}\}$, as defined in equation (42d). Recalling that this basis contains interleaved even harmonics of sines and cosines [see equations (17a) and (18a)]:

[TABLE]

we easily deduce that even elements $\breve{\beta}_{2p}=\sin(2p\theta)$ are antisymmetric functions whereas odd elements $\breve{\beta}_{2p-1}=\cos([2p-2]\theta)$ are symmetric functions. The symmetry of the operator $\mathscr{D}$ ensures us that this operator cannot couple symmetric and antisymmetric functions. Hence we only compute the coupling among (symmetric) cosines $\mathbfsf{D}^{c}$ and among (antisymmetric) sines $\mathbfsf{D}^{s}$ as products of elementary matrices $\mathcal{E}$ of size $N^{\prime}\geq N+1$:

[TABLE]

Finally, we obtain the expression of $\mathbfsf{D}$, the discretization of $\mathscr{D}=\sin^{2}\theta$ on the basis $\{\breve{\beta}_{n}\}$, by interweaving adequately the matrices $\mathbfsf{D}^{c}$ and $\mathbfsf{D}^{s}$:

[TABLE]

We indicate below the construction by matrix products for matrices of operators $\mathscr{A}=\sin^{6}\partial^{2}_{\theta\theta}$, $\mathscr{B}=\cos\theta\sin^{5}\theta\partial_{\theta}$, and $\mathscr{C}=\cos^{2}\theta$. These parity conserving operators are all treated in a similar fashion as $\mathscr{D}$: as justified above, we only compute coupling among cosines and sines:

[TABLE]

The final interweaving to obtain $\mathbfsf{A}$ and $\mathbfsf{B}$ is identical as in the case of $\mathbfsf{D}$, given in equations (103). The matrix $\mathbfsf{C}$ is directly obtained from the identity $\cos^{2}\theta=1-\sin^{2}\theta$. Thus, denoting $\mathbfsf{1}$ the identity matrix :

[TABLE]

The case of parity-flipping operators

Finally, we illustrate the treatment of a parity breaking term by discretizing the operator $\mathscr{E}=\cos\theta\sin\theta$. The parity-flipping nature of $\mathscr{E}$ guarantees that, in contrast with parity-conserving operators, this operator does not couple elements of the basis with the same parity. Instead, the coupling occurs between sines and cosines, leading us to compute $\mathbfsf{E}^{sc}$ and $\mathbfsf{E}^{cs}$:

[TABLE]

Finally, the matrix $\mathbfsf{E}$ defined in equation (42e) is obtained by the following interweaving, that differs from the parity-conserving case:

[TABLE]

XII.2 Discretization along a bounded direction: Chebyshev polynomials and the quasi-inverse technique

Central to the quasi-inverse technique is the three-term relation for the antiderivative of Chebyshev polynomials, first emphasized in the context of numerical resolution of ODEs by Clenshaw [23]:

[TABLE]

Upon integration, an arbitrary constant $k$ naturally appears. This constant will be set by enforcing boundary conditions: the consequence is that, for now, the top row of the matrix representing integration can be left equal to zero.

[TABLE]

Another handy three-term relationship corresponds to the multiplication by $x$:

[TABLE]

which has the discretized matrix:

[TABLE]

Hence the matrices $\mathbfsf{Q}$ defined in equations (76a,b,c) are again obtained by truncated matrix products of large matrices:

[TABLE]