Spectral Methods: Chebyshev Basis

Spectral Methods: Chebyshev Basis#

After the Fourier basis, then next most-useful basis is based on the Chebyshev polynomials, which can be expressed as

\[\begin{gather*} x = \cos\theta, \qquad T_n(x) = \cos(n\theta) = \cos(n\cos^{-1} x). \end{gather*}\]

A function expressed in this basis will thus have the form

\[\begin{gather*} f(x) = \sum_{n=0}^{N-1} a_n T_n(x) = \sum_{n=0}^{N-1} a_n \cos(n\theta) = F(\theta). \end{gather*}\]

Thus, expressing a function $f(x)$ for $x\in [-1, 1]$ in this basis is equivalent to considering a periodic function $F(\theta)$ with period $2\pi$ on the interval $\theta \in [0, \pi]$. These two representations are codified by the relationship:

\[\begin{gather*} x = \cos \theta, \qquad F(\theta) = f\underbrace{(\cos\theta)}_{x}. \end{gather*}\]

The relationship between Chebyshev and Fourier bases

If the Chebyshev basis is just a Fourier basis on $\theta$, where has the periodicity gone? The answer is that the periodic copies get mapped back into the interval $x\in [-1,1]$: the choice of $\theta \in [0, \pi]$ is because the mapping $\cos \theta$ is single-values / invertable on this interval.

This transformation implies that

\[\begin{gather*} F'(\theta) = - f'(x)\sin \theta, \qquad F''(\theta) = f''(x)\sin^2\theta - f'(x)\cos \theta,\\ f'(x) = -\frac{F'(\theta)}{\sin\theta}, \qquad f''(x) = \frac{F''(\theta) - \frac{\cos \theta}{\sin\theta}F'(\theta)}{\sin^2\theta}. \end{gather*}\]

Hence, if $f(x)$ is not singular, then $F'(\theta) = 0$ at the boundaries of the interval $[0, \pi]$ where $\sin \theta = 0$. Note however, that this does not imply any properties about the derivative $f'(x)$ at the boundaries.

The Chebyshev polynomials are orthogonal with respect to the weight $W(x) = 1/\sqrt{1-x^2} = 1/\sin\theta$:

\[\begin{gather*} \int_{-1}^{1} T_{n}(x)T_{m}(x)\frac{\d{x}}{\sqrt{1-x^2}} = \begin{cases} \pi & m=n=0\\ \tfrac{\pi}{2} & n=m \neq 0,\\ 0 & n \neq m. \end{cases} \end{gather*}\]

Integration and Differentiation#

Having a spectral representation affords us analytic integration and differentiation formulae

\[\begin{align*} T'_{n} &= nU_{n-1} = \sum_{j=0}^{n-1} c_j T_{n}, \qquad c_{0}, c_{1}, c_{2}, \cdots = \begin{cases} 0, 2, 0, 2, 2, 0, \cdots, & n\text{ even},\\ 1, 0, 2, 0, 2, 0, \cdots, & n\text{ odd},\\ \end{cases}\\ T''_{n} = nU_{n-1}' = \int T_{n} \equiv \int_{0}^{x} T_{n} \d{x} &= \frac{1}{2}\left( \frac{T_{n+1}}{n+1} - \frac{T_{n-1}}{n-1} \right) + \frac{n\sin(n\pi/2)}{n^2-1} \qquad (n\geq 2),\\ \int T_0 &= x = \frac{x + T_{-1}}{2}, \qquad \int T_1 = \frac{x^2}{2} = \frac{T_2}{4} + \frac{T_0}{4}, \qquad \end{align*}\]

Thus: if a function is exactly representable in the basis, then the derivative is exact, and the integral is almost exact except for the highest component.

We can use these as follows:

\[\begin{align*} f(x) &= \sum_{n=0}^{N}a_{n}\cos(n\theta), \\ f'(x) &= \frac{1}{\sin \theta} \sum_{n=1}^{N} na_{n} \sin(n\theta),\\ f''(x) &= \frac{-f'(x)\cos\theta}{\sin^2 \theta} - \frac{1}{\sin^2 \theta} \sum_{n=1}^{N} n^2a_{n} \cos(n\theta),\\ \int_{0}^{x} f(x)\d{x} &= a_0 \cos \theta + a_1 (\cos 2\theta + \cos 0\theta) +\\ &+ \sin\theta\sum_{n=2}^{N} \frac{-n a_n}{n^2-1}\sin(n\theta) + \cos\theta \sum_{n=2}^{N} \frac{-a_n}{n^2-1}\cos(n\theta)\\ &+ \sum_{n=2}^{N} \frac{n a_n}{n^2-1}\sin\frac{n\pi}{2}. \end{align*}\]

Each of these contains a discrete sine or cosine series that can be computed using the DCT and/or DST. For integration, it might be simpler just to shift the two series using the explicit formula, but this needs testing.

Some Details

\[\begin{gather*} x = \cos\theta, \qquad \d{x} = -\sin \theta \d{\theta}, \qquad \theta' = -1/\sin\theta, \\ T_{n}'(x) = -n\sin(n\theta)\theta' = \frac{n\sin (n \theta)}{\sin\theta}. \end{gather*}\]

\[\begin{gather*} \int T_{n}(x) \equiv \int_{0}^{x}T_{n}(x)\d{x} = -\int_{\pi/2}^{\theta} \sin\theta \cos(n\theta)\d{\theta}\\ = \frac{-n\sin θ\sin(nθ) - \cos θ\cos(nθ) + n\sin(\pi n/2))}{n^2 - 1}. \end{gather*}\]

To see this, note

\[\begin{align*} \diff{\sin \theta \sin(n\theta)}{\theta} &= -n\sin \theta \cos(n\theta) - \cos \theta \sin(n\theta)\\ \diff{\cos \theta \cos(n\theta)}{\theta} &= n\cos \theta \sin(n\theta) + \sin \theta \cos(n\theta),\\ \diff{n\sin \theta \sin(n\theta) + \cos \theta \cos(n\theta)}{\theta} &= (1-n^2)\sin \theta \cos(n\theta), \end{align*}\]

after canceling the wrong term. To check the previous formula, note

\[\begin{gather*} \frac{1}{2}\left( \frac{T_{n+1}}{n+1} - \frac{T_{n-1}}{n-1} \right) = \frac{(n-1)\cos[(n+1)\theta] - (n+1)\cos[(n-1)\theta]}{2(n^2-1)} = \\ = \frac{(n-1)\cos(n\theta + \theta) - (n+1)\cos(n\theta - \theta)]}{2(n^2-1)} = \\ = \frac{(n-1)[\cos\theta\cos(n\theta) - \sin\theta\sin(n\theta)] - (n+1)[\cos\theta\cos(n\theta) + \sin\theta\sin(n\theta)]}{2(n^2-1)} = \\ = \frac{-n\sin\theta\sin(n\theta) - \cos\theta\cos(n\theta)}{n^2-1} = \\ \end{gather*}\]

These can be manipulated to derive expressions in the basis, but it is easier to work algebraically using recurrence relationships. To get an idea about how these work, first introduce the complementary Chebyshev polynomials of the second kind

\[\begin{gather*} U_{n-1}(x) = \frac{\sin(n\theta)}{\sin\theta}. \end{gather*}\]

These satisfy the following recurrence relations:

\[\begin{align*} T_{0}(x) &= 1, & T_{1}(x) &= x, & T_{n+1}(x) &= 2xT_{n}(x) - T_{n-1}(x),\\ U_{0}(x) &= 1, & U_{1}(x) &= 2x, & U_{n+1}(x) &= 2xU_{n}(x) - U_{n-1}(x). \end{align*}\]

The first few of these are

\[\begin{align*} T_0 &= 1, & U_0 &= 1 &&= T_0\\ T_1 &= x, & U_1 &= 2x &&= 2T_1\\ T_2 &= 2x^2-1, & U_2 &= 4x^2 - 1 &&= T_0 + 2T_2\\ T_3 &= 4x^3-3x, & U_3 &= 8x^3 - 4x &&= 2(T_1 + T_3)\\ T_4 &= 8x^4-8x^2+1, & U_4 &= 16x^4 - 12x^2 + 1 &&= T_0 + 2T_2 + 2T_4\\ T_5 &= 16x^5 - 20x^3 + 5x, & U_5 &= 32x^5 - 32x^3 + 6x &&= 2(T_1 + T_3 + T_5). \end{align*}\]

To convert from one to the other:

\[\begin{align*} T_{n} &= \frac{U_{n} - U_{n-2}}{2},\\ U_{n} &= \begin{cases} \displaystyle 2\sum_{\text{odd } j>0}^{n} T_{j} = 2\sum_{i=0}^{(n-1)/2} T_{2i+1} & \text{for odd $n$},\\ \displaystyle 2\sum_{\text{even } j\geq 0}^{n} T_{j} - 1 = 2\sum_{i=0}^{n/2} T_{2i} - 1& \text{for even $n$}. \end{cases} \end{align*}\]

\[\begin{gather*} \int T_1 = \frac{T_{2} + T_{0}}{4} \int T_0 = \frac{T_{1}}{2} + \frac{T_{-1}}{2} = T_1 \end{gather*}\]

\[\begin{gather*} T'_{n}(x) = nU_{n-1}(x), \qquad U'_{n}(x) = \int T_{n}(x) = (1-n^2) \end{gather*}\]

functions

d2Tn_dx2 = dTn_dx.diff(th)/x.diff(th)
d2m = d2Tn_dx2.limit(th, sympy.pi).factor().simplify()
d2p = d2Tn_dx2.limit(th, 0).factor().simplify()
display(Latex(f"$T''_n(-1) = {sympy.latex(d2m)}$"))
display(Latex(f"$T''_n(+1) = {sympy.latex(d2p)}$"))

Pseudo-Spectral Chebyshev Basis#

The pseudo-spectral approach pairs a spectral basis $\phi_{m}(x)$ with an appropriate set of abscissa $x_{i}$. To be effective, one should be able to associate the grid with quadrature weights $w_i$ such that

\[\begin{gather*} \sum_{i}w_i f(x_i) = \int f(x) W(x)\d{x} \end{gather*}\]

is exact for polynomials $f(x)$ of degree $N$ or less (with $N$ abscissa).

There are two sets of quadrature points and weights based on the Chebyshev polynomials. Both are exact for polynomials of degree up to and including $x^{2N-1}$ (i.e., with $2N$ coefficients):

\[\begin{gather*} P(x) = a_0 + a_1x + a_2x^2 + \dots a_{2N-1} x^{2N-1}, \qquad \int_{-1}^{1} \frac{P(x)}{\sqrt{1-x^2}}\d{x} = \sum_{j} w_j P(x_j). \end{gather*}\]

We also include here the Cardinal Functions (DVR) $C_j(x)$ which are zero at all abscissa exact one

\[\begin{gather*} C_{j}(x_i) = \delta_{ij}. \end{gather*}\]

Chebyshev-Gauss (Interior Points / Roots)#

The $N$ abscissa are the roots of $T_{N}(x)$, which lie in the interior:

\[\begin{gather*} x_j = \cos \left.\left(\pi\frac{j + \tfrac{1}{2}}{N}\right)\right|_{j=0}^{N-1}, \qquad w_j = \frac{\pi}{N},\qquad C_{j}(x) = \frac{(x-x_j)T_{N}(x)}{T'_{N}(x_j)}. \end{gather*}\]

The $N$ weights and $N$ abscissa are free, allowing us to fit these $2N$ independent parameters. The relevant DCTs and DSTs are:

\[\begin{gather*} f_j = \sum_{k=0}^{N-1} a_{k} \cos(k\theta_j) = \sum_{k=0}^{N-1}b_{k}\sin\bigl((k+1)\theta_j\bigr),\\ \begin{alignedat}{2} [f_0, f_1, \cdots, f_{N-1}]& \overset{\text{DCT-III = iDCT-II}}{\underset{\text{DCT-II = iDCT-III}}{\rightleftharpoons}} &&[\tfrac{1}{2}a_{0}, \cdots, \tfrac{1}{2}a_{N-1}],\\ [f_0, f_1, \cdots, f_{N-1}]& \overset{\text{iDCT-III}}{\underset{\text{DST-III}}{\rightleftharpoons}} &&[\tfrac{1}{2}b_{0}, \cdots, \tfrac{1}{2}b_{N-2}, b_{N-1}]. \end{alignedat} \end{gather*}\]

Chebyshev-Gauss-Lobatto (Extrema + Endpoints)#

The $N+1$ abscissa are the extrema of $T_{N}(x)$ plus the endpoints:

\[\begin{gather*} x_j = \cos \left.\left(\pi\frac{j}{N}\right)\right|_{j=0}^{N}, \quad w_j = c_j\frac{\pi}{N},\\ C_{j}(x) = c_j\frac{(-1)^{j+1}(1-x^2)}{N^2}\frac{T'_{N}(x)}{x-x_j},\qquad c_j = \begin{cases} \frac{1}{2} & j\in \{0, N\},\\ 1 & \text{otherwise}. \end{cases} \end{gather*}\]

Since we fix the endpoints, we have $N+1$ weights and $N-1$ free abscissa, giving $2N$ free parameters. The relevant DCTs and DSTs are:

\[\begin{gather*} f_j = \sum_{k=0}^{N} a_{k} \cos\bigl(k\theta_j\bigr) = \left.\sum_{k=0}^{N-1} b_{k} \sin\bigl((k+1)\theta_{j}\bigr)\right|_{j=1}^{N},\\ \begin{alignedat}{2} [f_0, f_1, \cdots, f_{N}]& \overset{\text{iDCT-I}}{\underset{\text{DCT-I}}{\rightleftharpoons}} &&[a_{0}, \tfrac{1}{2}a_{1}\cdots, \tfrac{1}{2}a_{N-1}, a_{N}],\\ [f_1, f_2, \cdots, f_{N}]& \overset{\text{iDCT-III}}{\underset{\text{DST-III}}{\rightleftharpoons}} &&[\tfrac{1}{2}b_{0}, \cdots, \tfrac{1}{2}b_{N-1}]. \end{alignedat} \end{gather*}\]

These are implemented in the gpe.chebyshev.Chebyshev with the slight modification that the order of the points is changed so that the abscissa increase rather than decrease.

../_images/b3e3d9964c3fd8f30edaf4034efb6be5ce89277c9396662153446fc13526b9c6.png

Derivatives and Integration#

https://upload.wikimedia.org/wikipedia/commons/a/ae/DCT-symmetries.svg — Fig. 1 Periodic Extensions of the DCT#

DCT-I is even periodic about the endpoints, and hence suitable for the Gauss-Lobatto abscissa which have $x_{0} = -1$ and $x_{N} = 1$. DCT-2 is even about the external midpoints, and hence suitable for the Gauss abscissa.

https://upload.wikimedia.org/wikipedia/commons/3/31/DST-symmetries.svg — Fig. 2 Periodic Extensions of the DST#

DST-I is odd about the endpoints, and hence suitable for the Gauss-Lobatto abscissa which have $x_{0} = -1$ and $x_{N} = 1$. DST-2 is odd about the external midpoints, and hence suitable for the Gauss abscissa.

../_images/a0fd14298777f3833a3a9af97eace0f62888efba1e76370d5ba173ea546478d4.png

Discrete Cosine/Sine Transforms (DCT/DST)

As with the FFT, one can use the DCT and DST to compute derivatives. There are a variety of options differing in how the interval $\theta \in [0, \pi]$ is extended to a periodic function over $[0, 2\pi)$. DCT-I is periodic about the endpoints, and hence suitable for the Gauss-Lobatto abscissa which have; DCT-2 is periodic about the points outside of the see Periodic Extensions of the DCT. Here we discuss the appropriate formula for the DCT (see scipy.fft.dct()).

The strategy is:

Use the DCT to change bases, finding the coefficients $a_{k}$ such that

\[\begin{gather*} f_{j} = f(x_{j}) = \sum_{k}a_{k} T_{k}(x_{j}) = \sum_{k}a_{k} \cos (k \cos^{-1} x_{j}) = \sum_{k}a_{k} \cos (k \theta_{j}). \end{gather*}\]
Compute the derivative:

\[\begin{align*} f'(x_{j}) &= \sum_{k}a_{k} T_{k}'(x_{j}) = \sum_{k}a_{k} (-k)\sin(k \theta_{j})\theta'_{j},\\ &= \sum_{k}a_{k} k\frac{\sin(k \theta_{j})}{\sin \theta_{j}},\\ f''(x_{j}) &= \sum_{k}a_{k} T_{k}''(x_{j}) = \sum_{k}a_{k} (-k)\sin(k \theta_{j})\theta'_{j},\\ &= \sum_{k}a_{k}\left( -k^2\frac{\cos(k \theta_{j})}{\sin^2 \theta_{j}} + k\frac{\sin(k \theta_{j})\cos \theta_j}{\sin^3\theta_{j}} \right). \end{align*}\]
Use the DST and DCT to invert the transformation:

\[\begin{gather*} f_{j} = f(x_j) = \sum_{k}b_{k}\sin(k \theta_{j}) \end{gather*}\]

To do this in the code, we must choose the appropriate version of the DCT and DST, which differ in terms of the precise abscissa $\theta_{j}$ used, offsets in the indices, and factors for the coefficients. We first summarize the relevant transformations.

DCT-I:

\[\begin{align*} \theta_{k} &= \left.\frac{k\pi}{N}\right|_{k=0}^{N},\\ [f_0, f_1, \cdots, f_{N}] & \mathrel{\substack{\text{iDCT-I}\\\rightleftharpoons\\\text{DCT-I}}} [a_0, \tfrac{1}{2}a_{1}, \cdots, \tfrac{1}{2}a_{N-1}, a_{N}]. \end{align*}\]
DCT-II and DCT-III:

\[\begin{align*} \theta_{k} &= \left.\frac{(k+\tfrac{1}{2})\pi}{N}\right|_{k=0}^{N-1},\\ [f_0, f_1, \cdots, f_{N-1}] & \mathrel{\substack{\text{DCT-III = iDCT-II}\\\rightleftharpoons\\\text{DCT-II == iDCT-III}}} [\tfrac{1}{2}a_{0}, \cdots, \tfrac{1}{2}a_{N-1}]. \end{align*}\]
DST-I:

\[\begin{align*} \theta_{k} &= \left.\frac{k\pi}{N+1}\right|_{k=1}^{N},\\ [f_1, f_2, \cdots, f_{N}] & \mathrel{\substack{\text{DST-I}\\\rightleftharpoons\\\text{iDST-I}}} [\tfrac{1}{2}a_{1}, \cdots, \tfrac{1}{2}a_{N}]. \end{align*}\]
DST-III:

\[\begin{align*} \theta_{k} &= \left.\frac{k\pi}{N+1}\right|_{k=1}^{N},\\ [f_1, f_2, \cdots, f_{N}] & \mathrel{\substack{\text{DST-I}\\\rightleftharpoons\\\text{iDST-I}}} [\tfrac{1}{2}a_{1}, \cdots, \tfrac{1}{2}a_{N}]. \end{align*}\]

Type I:

From scipy.fft.dct(), changing only $N \rightarrow N+1$, $y_{k} \rightarrow f_{k}$, and $x_{n} \rightarrow \tilde{f}_{n}$ for notational convenience:

\[\begin{gather*} \underbrace{f_k}_{f(x_k)} = \underbrace{\tilde{f}_0}_{a_{0}}\underbrace{1}_{T_0(x_k)} + \underbrace{(-1)^{k}}_{T_{N}(x_k)} \underbrace{\tilde{f}_{N}}_{a_{N}} + \sum_{n=1}^{N-1} \underbrace{2\tilde{f}_{n}}_{a_{n}} \underbrace{\cos\Bigl(n\overbrace{\frac{\pi k}{N}}^{\theta_{k}}\Bigr)}_{T_{n}(x_{k})}. \end{gather*}\]

This matches the Gauss-Lobatto abscissa as shown with

\[\begin{gather*} \underbrace{\theta_{j}}_{\cos^{-1}(x_{k})} = \left.\frac{\pi j}{N}\right|_{j=0}^{N}, \qquad k \in \{0, 1, \dots, N\}. \end{gather*}\]

The DCT-I as written thus effects the inverse transformation $a_{n} \rightarrow f_{k}$, and the inverse effects the forward transform from $f_{k} \rightarrow a_{n}$ with appropriate scaling by factors of 2.

Similarly, the derivatives can be computed using the DST-I. Here is what we need to compute derivatives with respect to $\theta$, then we use the chain rule to compute $f'(x)$:

\[\begin{gather*} \underbrace{\diff{f_k}{\theta}}_{f'(\theta_k)} = - N\tilde{f}_{N}\sin(N\theta_{k}) - \sum_{n=1}^{N-1} 2n\tilde{f}_{n} \sin(n\theta_{k}). \end{gather*}\]

I.e., each coefficient $\tilde{f}_{n}$ is simply multiplied by $-n$, and now we have a sine series. To convert this back to position space, we use the DST-I, which computes

\[\begin{gather*} \underbrace{f_{k}}_{f(x_{k+1})} = \sum_{n=0}^{N-1} 2\tilde{f}_n \sin\Bigl((n+1)\overbrace{\frac{\pi(k+1)}{N+1}}^{\theta_{k+1}}\Bigr). \end{gather*}\]

Notice that we need to shift the coefficients. Finally, we need to use the chain rule to compute the derivative:

\[\begin{gather*} \underbrace{f'_{k}}_{f'(x_k)} = \frac{-1}{\sin\theta_k}\diff{f_k}{\theta_k} = N\tilde{f}_{N}\frac{\sin(N\theta_{k})}{\sin(\theta_k)} + \sum_{n=1}^{N-1} 2n\tilde{f}_{n} \frac{\sin(n\theta_{k})}{\sin(\theta_k)}. \end{gather*}\]

This is fine everywhere except at the endpoints where we must use L’Hôpital’s rule to obtain

\[\begin{gather*} \lim_{\theta \rightarrow 0}\frac{\sin(n\theta)}{\sin(\theta)} = n, \qquad \lim_{\theta \rightarrow \pi}\frac{\sin(n\theta)}{\sin(\theta)} = (-1)^{n-1}n,\\ f'(x_{0}) = f'(1) =\sum_{n=0}^{N-1} 2n^2\tilde{f}_{n} + N^2\tilde{f}_{N},\\ f'(x_{N}) = f'(-1) = \sum_{n=0}^{N-1} 2(-1)^{n-1} n^2\tilde{f}_{n} + (-1)^{N-1}N^2\tilde{f}_{N}. \end{gather*}\]

The second derivative has the following properties. Let $a_n = 2\tilde{f}_{n}$ except $a_{N} = \tilde{f}_{N}$. Then

\[\begin{gather*} f''(x) = \frac{f_{,\theta\theta}}{\sin^2\theta} - \frac{f_{,\theta}\cos\theta}{\sin^3\theta},\\ f'(\theta)= -\sum_{n=1}^{N} n a_n\sin(n\theta),\qquad f''(\theta)= -\sum_{n=1}^{N} n^2 a_n\cos(n\theta),\\ f''(x) = -\sum_{n=1}^{N} \frac{n^2\cos(n\theta)\sin\theta + n\sin(n\theta)\cos\theta}{\sin^3\theta}a_{n}. \end{gather*}\]

\[\begin{gather*} \lim_{\theta\rightarrow 0} \frac{n^2\cos(n\theta)\sin\theta - n\sin(n\theta)\cos\theta}{\sin^3\theta} = \frac{n^2(1-n^2)}{3},\\ \lim_{\theta\rightarrow \pi} \frac{n^2\cos(n\theta)\sin\theta - n\sin(n\theta)\cos\theta}{\sin^3\theta} = (-1)^{n}\frac{n^2(1-n^2)}{3}. \end{gather*}\]

Type III

\[\begin{gather*} \underbrace{f_k}_{f(x_{k})} = \underbrace{\tilde{f}_0}_{a_0} \underbrace{1}_{T_{0}(x_{k})} + \sum_{n=1}^{N-1} \underbrace{2\tilde{f}_n}_{a_{n}} \underbrace{ \cos\Bigl(n\overbrace{\frac{\pi(2k+1)}{2N}}^{\theta_{k}}\Bigr) }_{T_{n}(x_{k})} \end{gather*}\]

The DCT-II as written thus effects the inverse transformation $a_{n} \rightarrow f_{k}$, and the inverse transform (the DCT-II) effects the forward transform from $f_{k} \rightarrow a_{n}$ with appropriate scaling by factors of 2.

To compute derivatives, we note that

\[\begin{gather*} f'(\theta_{m}) = -2\sum_{n=0}^{N-1} \tilde{f}_{n} k_{n}\sin(k_{n}\theta_{m}),\\ f''(\theta_{m}) = -2\sum_{n=0}^{N-1} \tilde{f}_{n} k_{n}^2\cos(k_{n}\theta_{m}). \end{gather*}\]

These can be computed using the DST-III with coefficients $-\tilde{f}_{n}k_{n}$ shifted by one.

Similarly, the derivatives can be computed using the DST-II:

\[\begin{gather*} y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(2n+1)}{2N}\right) \end{gather*}\]

Similarly, the derivatives can be computed using the DST-III:

\[\begin{gather*} y_k = (-1)^k x_{N-1} + 2 \sum_{n=0}^{N-2} x_n \sin\left( \frac{\pi(2k+1)(n+1)}{2N}\right) \end{gather*}\]

Similarly, the derivatives can be computed using the DST-IV:

\[\begin{gather*} y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(2k+1)(2n+1)}{4N}\right) \end{gather*}\]

Here we demonstrate these formula:

import scipy.fft
sp = scipy

# Gauss quadrature at roots (interior points)

N = 12
n1 = k1 = np.arange(N)
th1 = np.pi * (2*n1 + 1) / 2 / N
x1 = np.cos(th1)

def get_a1(f):
    """Return the expansion coefficients."""
    ft = sp.fft.dct(f, type=2, norm='forward')
    a = 2*ft
    a[0] = ft[0]
    return a

def get_f1(a):
    """Return the function in position space."""
    ft = a/2
    ft[0] = a[0]
    return sp.fft.idct(ft, type=2, norm='forward')

def diff_1(f, d=1, th=th1, k=k1):
    ft = sp.fft.dct(f, type=2, norm='forward')
    s = np.sin(th)

    # Shift the coefficients and pad with zero.
    df_dth = sp.fft.dst(np.concatenate([-(ft*k)[1:], [0]]), type=3)
    if d == 1:
        df_dx = -df_dth / s
        return df_dx 
    elif d >= 2:
        d2f_dth2 = sp.fft.idct(-ft*k**2, type=2, norm='forward')
        d2f_dx2 = (d2f_dth2 - df_dth / np.tan(th)) / s**2
        if d == 2:
            return d2f_dx2
        return diff(d2f_dx2, d=d-2)

n2 = k2 = np.arange(N+1)
th2 = np.pi * n2 / N
x2 = np.cos(th2)

def get_a2(f):
    """Return the expansion coefficients."""
    ft = sp.fft.idct(f, type=1, norm='backward')
    a = 2*ft
    a[0], a[-1] = ft[0], ft[-1]
    return a

def get_f2(a):
    """Return the function in position space."""
    ft = a/2
    ft[0], ft[-1] = a[0], a[-1]
    return sp.fft.dct(ft, type=1, norm='backward')

def diff_2(f, d=1, th=th2, k=k2):
    ft = sp.fft.idct(f, type=1, norm='backward')
    s = np.sin(th)

    # Shift the coefficients and pad with zero.
    df_dth = sp.fft.dst(np.concatenate([-(ft*k)[1:], [0]]), type=1)
    if d == 1:
        df_dx = -df_dth / s
        return df_dx 
    elif d >= 2:
        d2f_dth2 = sp.fft.dct(-ft*k**2, type=1, norm='backward')
        d2f_dx2 = (d2f_dth2 - df_dth / np.tan(th)) / s**2
        if d == 2:
            return d2f_dx2
        return diff(d2f_dx2, d=d-2)

def T(x, n):
    return np.cos(n*np.acos(x))

# Test
x_ = np.linspace(-1, 1, 300)
for (x, n, get_a, get_f, diff) in [
    (x1, n1, get_a1, get_f1, diff_1),
    (x2, n2, get_a2, get_f2, diff_2)]:
    f = x * np.exp(-5*x**2)
    df = (1 - 10*x**2)*np.exp(-5*x**2)
    ddf = 10*x*(-3 + 10*x**2)*np.exp(-5*x**2)
    f_ = x_ * np.exp(-5*x_**2)
    df_ = (1 - 10*x_**2)*np.exp(-5*x_**2)
    ddf_ = 10*x_*(-3 + 10*x_**2)*np.exp(-5*x_**2)
    a = get_a(f)
    f1 = get_f(a)
    df1 = diff(f, d=1)
    ddf1 = diff(f, d=2)
    f2 = T(x[:, None], n[None, :]) @ a
    f1_ = T(x_[:, None], n[None, :]) @ a

    fig, axs = plt.subplots(1, 3, figsize=(10, 3), constrained_layout=True)
    ax = axs[0]
    ax.plot(x, f)
    ax.plot(x, f1)
    ax.plot(x, f2)
    ax.plot(x_, f_)
    ax.plot(x_, f1_)
    ax1 = ax.twinx()
    ax1.plot(x_, f1_ - f_, '--k')
    ax.set(xlabel="$x$", ylabel="$f(x) = xe^{-5x^2}$")
    ax1.set(ylabel="Error (dashed)");

    axs[1].plot(x, df)
    axs[1].plot(x, df1)
    axs[1].plot(x_, df_)
    axs[1].set(title=f"max err = {abs(df-df1).max():g}")

    axs[2].plot(x, ddf)
    axs[2].plot(x, ddf1)
    axs[2].plot(x_, ddf_)
    axs[2].set(title=f"max err = {abs(ddf-ddf1).max():g}");

/tmp/ipykernel_4831/3423782851.py:64: RuntimeWarning: divide by zero encountered in divide
  df_dx = -df_dth / s
/tmp/ipykernel_4831/3423782851.py:68: RuntimeWarning: divide by zero encountered in divide
  d2f_dx2 = (d2f_dth2 - df_dth / np.tan(th)) / s**2

../_images/eb79854e8c46fa8a37916d68e34198ee79fa3aa7e3bec8dade0f07fe09aa6a85.png

../_images/ff712b808e51770816a6dfdcd1f3c50af44a0c7b4a06cdba2f929fd38209ea6b.png

import scipy.fft
sp = scipy

# Gauss-Lobatto quadrature at extrema and endpoints

N = 12
n = k = np.arange(N+1)
th = np.pi * n / N
x = np.cos(th)

def get_a2(f):
    """Return the expansion coefficients."""
    ft = sp.fft.idct(f, type=1, norm='backward')
    a = 2*ft
    a[0], a[-1] = ft[0], ft[-1]
    return a

def get_f2(a):
    """Return the function in position space."""
    ft = a/2
    ft[0], ft[-1] = a[0], a[-1]
    return sp.fft.dct(ft, type=1, norm='backward')

def T(x, n):
    return np.cos(n*np.acos(x))

# Test
x_ = np.linspace(-1, 1, 300)
f = x * np.exp(-5*x**2)
f_ = x_ * np.exp(-5*x_**2)
a = get_a2(f)
f1 = get_f2(a)
f2 = T(x[:, None], n[None, :]) @ a
f1_ = T(x_[:, None], n[None, :]) @ a

fig, ax = plt.subplots()
ax.plot(x, f)
ax.plot(x, f1)
ax.plot(x, f2)
ax.plot(x_, f_)
ax.plot(x_, f1_)
ax1 = ax.twinx()
ax1.plot(x_, f1_ - f_, '--k')
ax.set(xlabel="$x$", ylabel="$f(x) = xe^{-5x^2}$")
ax1.set(ylabel="Error (dashed)");

../_images/de4463ace27f6fb86fc82b9f1524badb96e7018a758303db97756f76d1cdbdeb.png

Here is a simple implementation:

import scipy.fft
sp = scipy

Nx = 32
n = m = np.arange(Nx)
th = np.pi * (2*m + 1)/2/Nx
k = n
x = np.cos(th)

def diff(f, d=1, k=k, th=th):
    ft = sp.fft.dct(f, type=2, norm='forward')
    s = np.sin(th)

    # Shift the coefficients and pad with zero.
    df_dth = sp.fft.dst(np.concatenate([-(ft*k)[1:], [0]]), type=3)
    if d == 1:
        df_dx = -df_dth / s
        return df_dx 
    elif d == 2:
        d2f_dth2 = sp.fft.idct(-ft*k**2, type=2, norm='forward')
        d2f_dx2 = (d2f_dth2 - df_dth / np.tan(th)) / s**2
        return d2f_dx2
    return diff(df_dx, d=d-1)
    
f = np.exp(-5*x**2)
df = -10*x*f
ddf = (100 * x**2 - 10) * f

plt.plot(x, df, '-C0', alpha=0.5)
plt.plot(x, diff(f), '--C0')
plt.plot(x, ddf, '-C1', alpha=0.5)
plt.plot(x, diff(f, d=2), '--C1');

../_images/e701e44436aff1b25e775dc18e043429b6812002e03dee9066dda88f9926d46d.png

Cardinal Functions (DVR)#

The cardinal functions $C_j(x)$ are rearrangements of the basis functions such that

\[\begin{gather*} C_j(x_i) = \delta_{ij}. \end{gather*}\]

One can think of these as the delta-functions for the basis. To expand a function in the basis using these, one simply needs to evaluate the function at the abscissa:

\[\begin{gather*} f(x) \approx \sum_{j} f(x_j)C_j(x). \end{gather*}\]

Two sets of cardinal functions are given for the Chebyshev basis in [Boyd, 1989]. To evaluate these, we must be able to accurately compute

\[\begin{gather*} \frac{T_{N}'(x)}{x-x_j} = \frac{n\sin(n\theta)}{(\cos\theta - \cos\theta_j)\sin\theta}. \end{gather*}\]

Gauss-Lobatto Grid (Extrema)#

\[\begin{gather*} \left.x_j = \cos\frac{\pi j}{N}\right|_{j=0}^{N},\\ C_{j}(x) = \frac{(-1)^{j+1}(1-x^2)}{N^2}\frac{T'_{N}(x)}{x-x_j}\begin{cases} \frac{1}{2} & j\in \{0, N\},\\ 1 & \text{otherwise}. \end{cases} \end{gather*}\]

Roots Grid (Interior Points)#

\[\begin{gather*} \left.x_j = \cos\pi \frac{2j - 1}{2N}\right|_{j=1}^{N},\\ C_{j}(x) = \frac{(x-x_j)T_{N}(x)}{T'_{N}(x_j)}. \end{gather*}\]

def T(x, n, d=0):
    th = np.arccos(x)
    if d == 0:
        return np.cos(n*th)
    elif d == 1:
        return n*np.sin(n*th)/np.sin(th)
        
class Base:
    def __init__(self, **kw):
        for key in kw:
            if not hasattr(self, key):
                raise ValueError(f"Unknown {key=}")
            setattr(self, key, kw[key])
        self.init()
        
    def init(self):
        pass
    
class Basis1(Base):
    """Gauss-Lobatto Grid of extrema and endpoints."""
    N = 32
    def init(self):
        super().init()
        
    def get_x(self, j):
        return np.cos(np.pi * j / self.N)
        
    def get_C(self, x, j):
        N = self.N
        xj = self.get_x(j=j)
        C = (-1)**(j+1) * (1-x**2) / (x-xj) * T(x, N, d=1) / N**2 
        if j == N or j == 0:
            C /= 2
        return C

b = Basis1(N=4)
x = np.linspace(-1, 1, 256)
fig, ax = plt.subplots()
for j in range(0, b.N+1):
    l, = ax.plot(x, b.get_C(x, j=j), label=f"{j=}")
    ax.axvline([b.get_x(j=j)], c=l.get_c(), ls=":")
ax.legend();

/tmp/ipykernel_4831/3661239913.py:31: RuntimeWarning: invalid value encountered in divide
  C = (-1)**(j+1) * (1-x**2) / (x-xj) * T(x, N, d=1) / N**2
/tmp/ipykernel_4831/3661239913.py:6: RuntimeWarning: invalid value encountered in divide
  return n*np.sin(n*th)/np.sin(th)

../_images/1276e416cd657787ade14e80043f89db50631edd0c5f51bb501f5327ad198f7b.png

_EPS = np.finfo(float).eps

class Basis2(Base):
    """Interior points of Roots.
    
    Formulae from Appendix F of {cite}`Boyd:2001`.
    """
    N = 32
    
    def init(self):
        super().init()
        self.i = 1+np.arange(self.N)
        self.x = self.get_x(self.i)
        
    def get_x(self, j):
        return np.cos(np.pi * (2*j-1) / 2 / self.N)
        
    def get_C(self, x, j, d=0):
        N = self.N
        xj = self.get_x(j=j)
        t = np.arccos(x)
        tj = np.arccos(xj)
        return np.cos(N*t) * np.sin(tj) / N / np.sin(N*tj) / (np.cos(t) - np.cos(tj))

    def get_D1(self):
        """Return the derivative matrix."""
        i = self.i
        xi = self.x
        i, j = i[:, None], i[None, :]
        xi, xj = xi[:, None], xi[None, :]
        D1 = (-1)*(i+j)*np.sqrt((1-xj**2)/(1-xi**2))/(xi-xj+_EPS)
        D1 -= np.diag(np.diag(D1))
        xj = np.ravel(xj)
        D1 += np.diag(0.5*xj/(1-xj**2))
        return D1

    def get_D2(self):
        """Return the second derivative matrix."""
        N = self.N
        i = self.i
        xi = self.x
        i, j = i[:, None], i[None, :]
        xi, xj = xi[:, None], xi[None, :]
        D2 = self.get_D1() * (xi/(1-xi**2) - 2 / (xi-xj+_EPS))
        D2 -= np.diag(np.diag(D2))
        xj = np.ravel(xj)
        D2 += np.diag(xj**2/(1-xj**2)**2 - (N**2-1)/3/(1-xj**2))
        return D2

b = Basis2(N=4)
x = np.linspace(-1, 1, 256)
fig, ax = plt.subplots()
for j in range(1, b.N+1):
    l, = ax.plot(x, b.get_C(x, j=j), label=f"{j=}")
    ax.axvline([b.get_x(j=j)], c=l.get_c(), ls=":")
ax.legend();

../_images/8a95b7c84606d68a2106b4abbe5ca3a4f42befbab30ffafe7800a50e88b17a5b.png

\[T''_n(-1) = \frac{\left(-1\right)^{n} n^{2} \left(n^{2} - 1\right)}{3}\]

\[T''_n(+1) = \frac{n^{2} \left(n^{2} - 1\right)}{3}\]

A Chebyshev Basis with $f''(\pm 1)=0$.

One solution is linearly combine states to ensure this boundary condition. Note that

\[\begin{gather*} x = \cos\theta, \qquad T_n(x) = \cos(n\theta),\\ \theta'(x) = \frac{-1}{\sin\theta}, \qquad \theta''(x) = \frac{\theta'\cos\theta}{\sin^2\theta} = \frac{-\cos\theta}{\sin^3\theta}, \\ T_n'(x) = -n\sin(n\theta)\theta' = \frac{n\sin(n\theta)}{\sin\theta}, \\ T_n''(x) = n\frac{n\sin(n\theta)\cos\theta - n\cos(n\theta)\sin\theta}{\sin^3\theta}. \end{gather*}\]

Therefore:

\[\begin{align*} T_n(\pm 1) &= (\pm 1)^{n},\\ T_n'(\pm 1) &= (\pm 1)^{n}n^2,\\ T_n''(\pm 1) &= \frac{(\pm 1)^{n} n^2(n^2-1)}{3},\\ T_2''(\pm 1) &= 4,\\ T_3''(\pm 1) &= \pm 24. \end{align*}\]

Thus, $T_0$ and $T_1$ already have the correct boundaries. $T_2$ and $T_3$ don’t, so must be excluded, but can then be used to correct the basis functions. Here is a simple approach:

\[\begin{gather*} \varphi_{n}(x) = T_{n}(x) - \frac{n^2(n^2-1)}{12}T_2(x), \tag{even $n$}\\ \varphi_{n}(x) = T_{n}(x) - \frac{n^2(n^2-1)}{72}T_3(x). \tag{odd $n$} \end{gather*}\]

Unfortunately, these functions become highly degenerate because the corrections dominate for large $n$. Instead, we can do the following:

\[\begin{align*} \phi_{n}(x) &= T_{n}(x) - \frac{n^2(n^2-1)}{(n-2)^2((n-2)^2-1)} T_{n-2}(x)\\ &= T_{n}(x) - \frac{n^2(n+1)}{(n-2)^2(n-3)} T_{n-2}(x). \end{align*}\]

Thus, given

\[\begin{gather*} f(x) = \sum_{n=0}^{N-1} a_n T_n(x) = \sum_{n=0}^{N-1} b_n\phi_{n}(x)\\ = b_0 + b_1T_1(x) + \sum_{n=4}^{N-1} b_n T_{n}(x) - \sum_{n=4}^{N-1} \frac{b_n n^2(n+1)}{(n-2)^2(n-3)} T_{n-2}(x),\\ = b_0 + b_1T_1(x) + \sum_{n=4}^{N-1} b_n T_{n}(x) - \sum_{n=2}^{N-3} \frac{b_{n+2} (n+2)^2(n+3)}{n^2(n-1)} T_{n}(x),\\ \begin{multlined} = b_0 + b_1T_1(x) - 20 b_{4} T_{2}(x) - \frac{25}{3}b_{5} T_{3}(x)\\ + \sum_{n=4}^{N-3} \left(b_n - \frac{b_{n+2} (n+2)^2(n+3)}{n^2(n-1)}\right) T_{n}(x) + \sum_{n=N-2}^{N-1} b_n T_{n}(x). \end{multlined} \end{gather*}\]

This gives us the following recurrence

\[\begin{gather*} a_0 = b_0\\ a_1 = b_1\\ a_2 = -20b_4\\ a_3 = -\frac{25}{3}b_5\\ \left.a_n = b_n - b_{n+2}\frac{(n+2)^2(n+3)}{n^2(n-1)}\right|_{n=4}^{N-3}\\ a_{N-2} = b_{N-2}\\ a_{N-1} = b_{N-1}. \end{gather*}\]

\[\begin{gather*} b_{N-1} = a_{N-1},\\ b_{N-2} = a_{N-2},\\ \left.b_n = a_n + b_{n+2}\frac{(n+2)^2(n+3)}{n^2(n-1)}\right|_{n=4}^{N-3}\\ b_5 = -\frac{3}{25}a_3\\ b_4 = \frac{a_2}{20}\\ b_1 = a_1\\ b_0 = a_0\\ \end{gather*}\]

Nx = 32
n = m = np.arange(Nx)
th = np.pi * (2*m + 1) / 2/ Nx
k = n
x = np.cos(th)

f = np.sin(np.pi * x**2)
a = get_a1(f)
b = a.copy()
for n in range(Nx-3, 3, -1):
    b[n] = -a[n] + b[n+2]*(n+2)**2*(n+3)/n**2/(n-1)
b[2] = b[3] = 0

x = np.linspace(-1, 1, 500)
def T(x, n):
    return np.cos(n*np.acos(x))

def phi1(x, n):
    fact = n**2*(n**2 - 1)
    
    if n % 2 == 0:
        phi = T(x, n) - fact*T(x, 2)/12
    else:
        phi = T(x, n) - fact*T(x, 3)/72
    return phi / phi.max()

def phi2(x, n):
    fact = n**2*(n**2 - 1)/(n-2)**2/((n-2)**2 - 1)
    phi = T(x, n) - fact*T(x, n-2)
    return phi / abs(phi).max()

def d2(f):
    for n in range(2):
        f = np.gradient(f, x, edge_order=2)
    return f
    
for n in range(4, 10):
    l, = plt.plot(x, phi2(x, n), label=f"{n=}")
    d2_phi = d2(phi2(x, n))
    plt.plot(x, d2_phi/abs(d2_phi).max(), ':', c=l.get_c())
plt.legend();

../_images/48fd0394b2f691e30e4ceba403a0a51a6cd99092874fd5705eabc9a1a635f3d6.png

The first spectral basis we will consider is provided by the Chebyshev polynomials of the first kind $\{T_n(x)\}$ which are orthogonal on the interval $x\in [-1, 1]$ with weighting function $W(x) = 1/\sqrt{1-x^2}$:

\[\begin{gather*} T_{n}(\cos \theta) = \cos(n\theta), \qquad T_{n}(x) = \cos(n\cos^{-1}x),\\ \int_{-1}^{1} T_{m}(x)T_{n}(x)\frac{\d{x}}{\sqrt{1-x^2}} = \begin{cases} \pi & n = m = 0,\\ \frac{\pi}{2} & n = m \neq 0,\\ 0 & n \neq m. \end{cases} \end{gather*}\]

The Chebyshev polynomials of the second kind $U_{n}(x)$ are useful for computing derivatives:

\[\begin{gather*} U_n(\cos\theta) = \frac{\sin\bigl((n+1)\theta\bigr)}{\sin\theta},\\ T_n' = nU_{n-1}, \qquad U_n' = \frac{(n+1)T_{n+1} - xU_n}{x^2 - 1},\\ T_n'' = n\frac{nT_{n} - xU_{n-1}}{x^2 - 1} = n\frac{(n+1)T_{n} - U_n}{x^2 - 1}. \end{gather*}\]

If we express

\[\begin{gather*} f(x) = \sum_{n} f_nT_n(x) \end{gather*}\]

from scipy.integrate import quad

def T(x, n):
    return np.cos(n*np.acos(x))
    
quad(lambda x: T(x,2)*T(x,2)/np.sqrt(1-x**2), -1, 1)

(1.570796326794891, 1.7586718303874704e-08)

Spectral Methods: Chebyshev Basis

Contents

Spectral Methods: Chebyshev Basis#

Integration and Differentiation#

Pseudo-Spectral Chebyshev Basis#

Chebyshev-Gauss (Interior Points / Roots)#

Chebyshev-Gauss-Lobatto (Extrema + Endpoints)#

Refinement#

Derivatives and Integration#

Cardinal Functions (DVR)#

Gauss-Lobatto Grid (Extrema)#

Roots Grid (Interior Points)#