--- execution: timeout: 60 jupytext: cell_metadata_filter: -all formats: md:myst,ipynb main_language: python text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.7 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os, sys from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None from matplotlib.animation import FuncAnimation from gpe.contexts import FPS ``` :::{margin} The Chebyshev basis is just a Fourier basis on $\theta$ (with only even terms) disguised by a coordinate transformation from $\theta$ to $x = \cos^{-1}\theta$. ::: (sec:ChebyshevBasis)= # Spectral Methods: Chebyshev Basis After the Fourier basis, then next most-useful basis is based on the Chebyshev polynomials, which can be expressed as :::{margin} For practical purposes, we allow the original interval to be scaled and translated: \begin{gather*} x = L\cos\theta + x_L. \end{gather*} This gives simple corrections to the formulae. ::: \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*} :::{admonition} The relationship between Chebyshev and Fourier bases :class: dropdown 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. ::: :::{margin} This weight just comes from changing variables in the Fourier orthogonality relation \begin{gather*} \int_{0}^{\pi} \cos(n\theta)\cos(m\theta)\d{\theta}\\ = \int_{1}^{-1} T_{n}(x)T_{m}(x)\underbrace{\diff{\theta}{x}}_{-1/\sin\theta}\d{x}\\ = \int_{-1}^{1} T_{n}(x)T_{m}(x)\frac{\d{x}}{\sqrt{1-x^2}}. \end{gather*} ::: 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*} :::{margin} The integration formula works for $n=0$ if we take $T_{-1} = x$. For $n=1$, we must just drop the divergent term $\int T_1 = T_2/4$. ::: ```{code-cell} :tags: [margin, hide-cell] import sympy theta, x = sympy.symbols('theta, x') n = sympy.symbols('n', integer=True, positive=True) # Define based on recursive formula T = [1, x] U = [1, 2*x] for n in range(1, 10): T.append((2*x*T[n] - T[n-1]).simplify().expand()) U.append((2*x*U[n] - U[n-1]).simplify().expand()) print(",\\\\\n".join( [f"T_{{{n}}} &= {sympy.latex(T)}, & U_{{{n}}} &= {sympy.latex(U)}" for n, (T, U) in enumerate(zip(T, U))]) ) for n in range(2, 10): print(n, (sympy.integrate(T[n], (x, 0, x)) - (T[n+1]/(n+1) - T[n-1]/(n-1))/2).simplify()) print(sympy.sin(n*sympy.pi/2)*n/(n**2-1)) ``` ## 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. :::{admonition} 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 ::: ```{code-cell} :tags: [margin, hide-input] from IPython.display import Latex import sympy from sympy import var, cos, sin th = var('theta') n = var('n', integer=True) Tn = cos(n*th) Un = sin((n+1)*th)/sin(th) x = cos(th) dTn_dx = Tn.diff(th)/x.diff(th) assert dTn_dx == n*Un.subs([(n, n-1)]) ``` ``` 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 :::{margin} For orthogonal polynomials $p_{n}(x)$ ($n\in \{0, 1, \dots, N-1$), one can choose both $N$ weights $w_i$ and $N$ abscissa $x_i$ to make this exact for polynomials of degree $2N-1$. The strategy is to use the roots of $p_{N}(x_i) = 0$, then to compute the weights with \begin{gather*} w_{j} = \frac{\braket{p_{N-1}|p_{N-1}}}{p_{N-1}(x_j)p'_{N}(x_j)}, \tag{4.5.9}. \end{gather*} (From {cite}`PTVF:2007`.) ::: 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 [DCT][]s and [DST][]s 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 [DCT][]s and [DST][]s 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 {class}`gpe.chebyshev.Chebyshev` with the slight modification that the order of the points is changed so that the abscissa increase rather than decrease. ```{code-cell} :tags: [hide-input, margin] from gpe.chebyshev import Chebyshev c1 = Chebyshev(N=16) f1 = np.exp(-5*c1.x**2) c2 = Chebyshev(N=32) # Here we use the basis to extrapolate to the finer grid. Internally # c2.get_f() zero-pads (or truncates) the coefficients. f2 = c2.get_f(c1.get_a(f1)) fig, ax = plt.subplots(figsize=(3, 2)) ax.plot(c1.x, f1, label=f"$N_x={c1.N}$") ax.plot(c2.x, f2, label=f"$N_x={c2.N}$") ax.set(xlabel="$x$", ylabel="$e^{-5x^2}$") ax.legend(); ``` ### Refinement Refining in the Chebyshev basis is trivial: to start resolving more features, simply add higher-order coefficients, setting them to zero. :::::{admonition} Chebyshev Quadrature (Interior Points) The roots of $T_{N}(x_j) = 0$ are trivial: \begin{gather*} T_{N}(x_j) = \cos\bigl(N\theta_j\bigr) = 0, \qquad T'_{N}(x_j) = \frac{N\sin(N\theta_j)}{\sin\theta_j},\\ \theta_j = \frac{\pi (2j + 1)}{2N}, \qquad x_j = \cos\theta_j = \left.\cos\left(\pi\frac{2j + 1}{2N}\right)\right|_{j=0}^{N-1}. \end{gather*} The weights are thus \begin{gather*} w_j = \frac{\pi p_{N-1} \sin\theta_j} {2 N \cos\bigl((N-1)\theta_{j}\bigr)\sin(N\theta_j)}, \\ p_{0} = 2, \qquad p_{n>0} = 1. \end{gather*} :::{admonition} Numerical Check Here we check the accuracy of these quadrature schemes using \begin{gather*} \braket{x^n} = \int_{-1}^{1}\d{x}\; \frac{x^{n}}{\sqrt{1-x^2}} = \int_{0}^{\pi}\d{\theta} \cos^{n}\theta = \begin{cases} \frac{\pi}{2^{n/2-1}}\binom{n-1}{n/2} & n \text{ even},\\ 0 & n\text{ odd}. \end{cases} \end{gather*} *(I found this from .)* ``` import sympy th = sympy.var('theta') [sympy.integrate(sympy.cos(th)**(2*n), (th, 0, sympy.pi)) * 2**(2*n-1) / sympy.pi for n in range(1, 12)] ``` ::: ::::: ```{code-cell} import sympy th = sympy.var('theta') [sympy.integrate(sympy.cos(th)**(2*n), (th, 0, sympy.pi)) * 2**(2*n-1) / sympy.pi for n in range(1, 12)] ``` ```{code-cell} # Chebyshev Quadrature from functools import partial import scipy.integrate import scipy as sp def get_xw1(N): """Return the N Chevbyshev absicssa and weights (interior points).""" j = np.arange(N) th = np.pi * (2*j + 1)/2/N x = np.cos(th) w = np.pi / N return x, w def get_xw2(N): """Return the N+1 Chevbyshev-Lobatto absicssa and weights (extrema+endpoints).""" j = np.arange(N+1) th = np.pi * j / N x = np.cos(th) w = np.pi / N + 0*x w[0] /= 2 w[-1] /= 2 return x, w # Three ways to check # Sympy: This uses hypergeometric functions import sympy x = sympy.var('x') n = sympy.var('n', integer=True, positive=True) I_sympy = sympy.lambdify([n], sympy.integrate(x**n/sympy.sqrt(1-x**2), (x, -1, 1)), 'mpmath') # Exact formula def I_exact(n): if n == 0: return np.pi elif n %2 == 1: return 0 elif n % 2 == 0: return np.pi * sp.special.binom(n - 1, n//2) / 2**(n-1) # Adaptive quadrature def W(x): return 1/np.sqrt(1-x**2) def integrand(x, n): return x**n * W(x) def I_quad(n): res, err = sp.integrate.quad( partial(integrand, n=n), -1, 1, epsabs=1e-11, epsrel=1e-12) return res, err ``` ```{code-cell} for N in range(1, 10): ns = np.arange(2*N+1) cheb1 = [] cheb2 = [] for n in ns: x, w = get_xw1(N) cheb1.append((w*x**n).sum()) x, w = get_xw2(N) cheb2.append((w*x**n).sum()) exact = np.array([I_exact(n) for n in ns]) quad = [I_quad(n) for n in ns] sym = [I_sympy(n) for n in ns] assert np.allclose(exact[:-1], cheb1[:-1]) assert np.allclose(exact[:-1], cheb2[:-1]) assert not np.allclose(exact[-1], cheb1[1]) assert not np.allclose(exact[-1], cheb2[1]) fig, ax = plt.subplots(figsize=(5, 2)) ax.semilogy(ns, abs(exact - cheb1), label="Chebyshev-Gauss (interior)") ax.semilogy(ns, abs(exact - cheb1), label="Chebyshev-Lobatto (extrema + endpoints)") ax.set(xlabel="$n$", xticks=ns, ylabel="Absolute error", title=f"Error integrating $x^n$ with ${N=}$") ax.legend(); ``` ### Derivatives and Integration :::{figure} https://upload.wikimedia.org/wikipedia/commons/a/ae/DCT-symmetries.svg --- figclass: margin name: fig:DCT --- 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. ::: :::{figure} https://upload.wikimedia.org/wikipedia/commons/3/31/DST-symmetries.svg --- figclass: margin name: fig:DST --- 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. ::: ```{code-cell} :tags: [hide-input, margin] import timeit import scipy as sp rng = np.random.default_rng(seed=2) ns = {} ts = {} Nmax = 2**14 for p in [2, 3, 5]: ns[p] = p**np.arange(4, 1+int(np.log(Nmax)/np.log(p))) ns['prime'] = list(sympy.ntheory.generate.primerange(2**4, Nmax)) ns['prime'] = np.array(ns['prime'][::(len(ns['prime'])//20)]) #ns['2 prime'] = 2 * (ns['prime'] // 2) ns['4 prime'] = 4 * (ns['prime'] // 4) ts = {key: {} for key in ns} Nr = 20 # We do these in a random order so that the variations are more meaningful. # Otherise, all results from a particular run might be contaminated by a # background process. tasks = rng.permutation(np.array([ (key, Nx, r) for key in ns for Nx in ns[key] for r in range(Nr)], dtype=object)) for (key, Nx, r) in tasks: x = np.cos(np.linspace(-1, 1, Nx)) f = np.exp(-x**2) ts[key].setdefault(Nx, []).append( timeit.timeit('sp.fft.dct(f)', globals=globals(), number=1)) fig, ax = plt.subplots(figsize=(4, 4)) for lab in ns: n = sorted([key for key in ts[lab]]) t_min = np.array([np.min(ts[lab][key]) for key in n]) t_med = np.array([np.median(ts[lab][key]) for key in n]) l, = ax.plot(n, t_min*1e6, '-o', label=lab if isinstance(lab, str) else f"${lab}^n$") ax.plot(n, t_med*1e6, '--', c=l.get_c()) ax.set(xlabel="$N_x$", ylabel="time [μs]") plt.legend(); ``` :::{margin} This figure shows the cost of a DCT/iDCT pair on a function of length $N_x$, demonstrating that powers of small primes $N_x = p^n$ perform better than primes. Each curve has the minimum of 20 tests as a symbol with the median as a dashed line to allow you to estimate variability on your machine. ::: :::{admonition} 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 {ref}`fig:DCT`. Here we discuss the appropriate formula for the [DCT][] (see {func}`scipy.fft.dct`). The strategy is: 1. 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*} 2. 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*} 3. 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 {func}`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: ```{code-cell} 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}"); ``` ```{code-cell} 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)"); ``` Here is a simple implementation: ```{code-cell} 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'); ``` ## 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 {cite}`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*} ```{code-cell} 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(); ``` ```{code-cell} _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(); ``` ```{code-cell} :tags: [margin, hide-input] from IPython.display import Latex import sympy th = sympy.var('theta') n = sympy.var('n', integer=True) Tn = sympy.cos(n*th) x = sympy.cos(th) dTn_dx = Tn.diff(th)/x.diff(th) 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)}$")) ``` :::{admonition} 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*} ::: ```{code-cell} 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 ``` ```{code-cell} 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(); ``` [DCT]: [DST]: 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*} ```{code-cell} 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) ``` [Bessel function]: [Bohr radius]: [Jacobi elliptic functions]: [Jupyter Book with Sphinx]: [Jupyter Book]: [Jupyter]: "Jupyter" [Jupytext]: "Jupyter Notebooks as Markdown Documents, Julia, Python or R Scripts" [Laplace-Beltrami operator]: [Liouville's Theorem]: [Manim Community]: [Markdown]: [MyST Cheatsheet]: [MyST]: "MyST - Markedly Structured Text" [MySt-NB]: [Rutherford scattering]: [Sphinx]: [angular momentum operator]: [glue]: [hydrogenic atoms]: [Gaussian quadrature]: [Dirac delta function]: [compact support]: [hydrogen atom]: [reduced Bohr radius]: [generalized Laguerre polynomial]: [spherical harmonic]: [finite difference methods]: [Dirichlet boundary conditions]: [numerical integration]: [DVR]: [spectral methods]: [FEM]: [orthogonal polynomials]: [principal quantum number]: [azimuthal quantum number]: [magnetic quantum number]: [Numba]: [Whittaker-Shannon interpolation formula]: [Gauss-Radau rules]: [Gauss-Lobatto rules]: [Laguerre polynomials]: [generalized Laguerre polynomials]: [Lebesque-Stieltjes integral]: [L'Hôpital's rule]: