Chebyshev Transforms

A necessary side product of the conversion of the SWT code to use FFTW was a routine to carry out Chebyshev transforms. Had the FFT API of 2DECOMP&FFT been used, it would have been necessary to incorporate this into that library; as it stands, the domain decomposition API was chosen instead, so this is not necessary. However, it was considered appropriate to adapt and generalise this routine so as to make it suitable for use in or with 2DECOMP&FFT. This would broaden the capabilities of this library by reusing a side-product of this project.

Given a function $f(x)$ on the interval $-1 \le x \le 1$, and making a change of variable

$$x = \cos(\theta),$$ (12)

where $0 \le \theta \le \pi$, we define $\phi(\theta) = f(\cos(\theta))$. $\phi$ can be expanded in a cosine series

$$\phi(\theta) = \frac{1}{2} a_0 + \sum_{k=1}^\infty a_k \cos(k \theta)$$ (13)

The Chebyshev polynomials $T_k$ of the first kind are defined

$$T_k(x) = T_k(\cos(k \theta)) = T_k(k \arccos(x))$$ (14)

so

$$\phi(\theta) = \frac{1}{2} a_0 + \sum_{k=1}^\infty a_k T_k(x)$$ (15)

Note that the coefficients $a_k$ in 13 and 15 are the same. It is, therefore, possible to obtain them using a cosine transform, computable with the FFT algorithm. It is also possible, with some manipulation, to use a standard FFT. The SWT code originally used a complex-to-real FFT for this purpose, and Chebyshev transform routines in $x$, $y$ and $z$ (or rather the array indeces usually corresponding to those dimensions) were developed from this. However, as FFTW3 supports cosine transforms, an alternative set of routines was created taking advantage of these. Both sets have been made available to Dr. Ning Li, the developer of 2DECOMP&FFT.

Collocation methods based on Chebyshev polynomials usually use the Chebyshev-Gauss points located at the roots of $T_k$, and given by

$$x_k = \cos \left( \frac{2k+1}{2N} \pi \right), k=0,1..,N$$ (16)

or the Gauss-Lobatto points, located at the extrema of $T_k$ (in other words, at the zero points of $T_k\prime$ rather than $T_k$). These are given by

$$x_k = \cos \left( \frac{k}{N} \pi \right), k=0,1..,N$$ (17)

For applications such as SWT - which is designed to simulate flow through a channel of fixed depth - the Gauss-Lobatto points are preferred; the zeroth and Nth points are at 1 and -1, and do not depend, as the Chebyshev-Gauss points do, on $N$. Both grids are supported; at present, the default is to use the Gauss-Lobatto points, and an optional argument is used (if present) to switch to Chebyshev-Gauss.