--- 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 ``` # Spectral Methods Our general goal is to find an efficient basis for representing our functions: \begin{gather*} f(x) = \sum_{n} a_n \phi_n(x). \end{gather*} The idea is to then choose the coefficients $a_n$ so that $f(x)$ is as close as possible to the solution we desire -- e.g., solving a differential equation, or minimizing an energy functional. :::{admonition} Some Terminology The term **[spectral method][spectral methods]** applies to expansions such as those above where the functions $\phi_{n}(x)$ are generally analytic -- often polynomials. If the basis functions have [compact support][], then such an expansion this is referred to as a [finite-element method (FEM)][FEM]. The [FEM][] approach is often better for systems with discontinuities, irregular boundaries, and sharp features like shockwaves, but lacks the exponential convergence enjoyed by [spectral methods][] for analytic functions. The term **pseudo-spectral method**, also known as a [discrete variable representation (DVR)][DVR] methods, associates a quadrature grid with the basis functions, and reorganize the basis into a set of **cardinal functions** which vanish on all but one of the grid points. These are the representation of the [Dirac delta function][] in the basis, and greatly facilitate the expression of problems. ::: :::{admonition} Literature (How I learned about this.) :class: dropdown The Wikipedia entry about [spectral methods][] provides a good overview of the relationship between [spectral methods][], [finite element methods][FEM], [pseudo-spectral methods (DVR)][DVR] and other related methods. The following two textbooks provide a fairly comprehensive starting point for learning about spectral methods. * {cite}`Boyd:1989` -- "Chebyshev and Fourier Spectral Methods": This is an excellent resource, although the writing is quite idiosyncratic. Start with the introduction to get a very hands-on and practical introduction. The rest of the book is quite difficult to read, but as you learn more, keep coming back. There is a lot of practical wisdom here, and many references to real-world applications. * {cite}`Trefethen:2000` -- "Spectral Methods in Matlab": This is a little easier to read than {cite}`Boyd:1989`, with many concrete examples, but is not as thorough or deep. For applications in quantum mechanics, and my first exposure to the idea, see the following papers: * {cite}`Baye:1986si` -- "Generalized meshes for quantum mechanical problems": Demonstrates the essential features of [DVR][] bases for solving quantum mechanics problems. * {cite}`LCCMP:2002` -- "A general framework for discrete variable representation basis sets": This presents arguments about why DVR bases work so well in terms of phase-space coverage, and discusses many standard sets. * {cite}`LC:2002` -- "Bessel discrete variable representation bases": Go here if you need the formulae and theory for spherically symmetric problems in various dimensions (i.e., cylindrical or spherical geometries). * {cite}`Schneider:2004` -- "Discrete variable representation for singular Hamiltonians": I really enjoyed this presentation. It was the first time I understood the generalized relationship between quadrature and DVR bases, with a clever but simple generalization for dealing with singularities like the Coulomb potential. ::: ## Background / Prerequisites This section requires the following concepts: * **Orthonormal functions as basis vectors**: A set of functions $P_{n}(x)$ are said to be **orthonormal** with respect to an integration weight $w(x)$ iff \begin{gather*} \braket{P_{m}|P_{n}} \equiv \int w(x)\d{x}\; P_{m}(x)P_{n}(x) = \delta_{mn}. \end{gather*} * **(Gaussian) Quadrature**: An $N$-point quadrature rule is a set of $N$ abscissa $x_{n}$ and weights $w_{n}$ for $n \in \{0, 1, \dots, N-1\}$ such that \begin{gather*} \int w(x)\d{x}\; f(x) = \sum_{n=0}^{N-1} w_n f(x_n) \end{gather*} is exact for all $f(x) = P_{n}(x)$ with $n < N$. The term [Gaussian quadrature][] refers to a set of abscissa and points such that the integration formula is exact for all *polynomials* of degree $2N-1$. This is consistent with there being $2N$ degrees of freedom -- the weights and the abscissa -- and that there are $2N$ independent monomials $x^0$, $x^1$, ... $x^{2N-1}$. * **Fourier Transform and Projections**: You should also be familiar with the Fourier transform and projecting functions onto states of low momentum. * **Dirac Delta Functions**: The most general approach is in terms of $\delta(x)$. ## Paths to a [DVR][] ### TL;DR The best exposition I have seen of [DVR][] bases is from {cite}`Meyer:2011` which expresses the following theorem: :::::{prf:theorem} DVR Bases Consider a pseudo-spectral bases of $N$ orthonormal functions $P_n$, abscissa $x_a$, and quadrature weights $w_a$. Let \begin{gather*} [\mat{G}]_{an} = P_{n}(x_a), \qquad [\mat{w}]_{ab} = w_{a}\delta_{ab}, \qquad \mat{U} = \mat{G}^\dagger\mat{w}^{1/2},\\ \ket{F_n} = \sum_{n}\ket{P_m}U_{mn}. \end{gather*} The following 5 statements are *equivalent*: 1. $\mat{G}\mat{G}^\dagger = \mat{w}^{-1}$ is diagonal. 2. $\mat{U}^\dagger = \mat{U}^{-1}$ is unitary. 3. Discrete orthogonality: \begin{gather*} \sum_{a=0}^{N-1} w_a P_i^*(x_a)P_j(x_a) = \delta_{ij} \end{gather*} 4. Discrete completeness: \begin{gather*} \sum_{n=0}^{N-1} w_a P_n^*(x_a)P_n(x_b) = \delta_{ab} \end{gather*} 5. Discrete $\delta$-property (locality): \begin{gather*} F_n(x_a) = \frac{\delta_{na}}{\sqrt{w_n}}, \qquad \braket{F_m|F_n} = \delta_{mn}. \end{gather*} :::{prf:proof} :class: dropdown Much of this follows trivially when expressing everything in braket notation. To this end, let \begin{gather*} \ket{\mat{P}} = \begin{pmatrix} \ket{P_0}\\ \ket{P_1}\\ \vdots\\ \ket{P_{N-1}} \end{pmatrix}, \qquad \ket{\mat{x}} = \begin{pmatrix} \ket{x_0}\\ \ket{x_1}\\ \vdots\\ \ket{x_{N-1}} \end{pmatrix}, \end{gather*} In terms of these, we define the following \begin{gather*} \mat{G} = \braket{\mat{x}|\mat{P}}, \qquad \mat{U} = \braket{\mat{P}|\mat{x}}\mat{w}^{1/2}, \qquad \ket{\mat{F}} = \ket{\mat{P}}\mat{U}, \end{gather*} where $\mat{w}$ is real, diagonal, and invertable. The orthonormality of $\ket{P_n}$ means \begin{gather*} \braket{P_m|P_n} = \delta_{mn}, \qquad \braket{\mat{P}|\mat{P}} = \mat{1}. \end{gather*} Then \begin{gather*} \mat{U}^\dagger\mat{U} = \mat{w}^{1/2}\mat{G}\mat{G}^\dagger\mat{w}^{1/2}. \end{gather*} So if $\mat{G}\mat{G}^\dagger = \mat{w}^{-1}$ is diagonal, then $\mat{U}$ is unitary ($1. \Rightarrow 2.$), and if $\mat{U}$ is unitary, then we can multiply both sides by $\mat{w}^{-1/2}$ to show ($2. \Rightarrow 1.$). We now just manipulation the expressions: \begin{gather*} \mat{1} = \mat{G}\mat{G}^\dagger\mat{w} = \braket{\mat{x}|\mat{P}}\braket{\mat{P}|\mat{x}}\mat{w} \end{gather*} Property 2. implies that the DVR basis $\ket{F_n}$ is orthonormal and local: \begin{gather*} \braket{\mat{x}|\mat{F}} = \braket{\mat{x}|\mat{P}}\mat{U} = \mat{w}^{-1/2}\mat{U}^\dagger\mat{U} = \mat{w}^{-1/2},\\ \braket{\mat{F}|\mat{F}} = \mat{U}^\dagger\braket{\mat{P}|\mat{P}}\mat{U} = \mat{U}^\dagger\mat{U} = \mat{1}. \end{gather*} Hence $2. \Leftrightarrow 5.$ The remaining properties 3. and 4. follow from the fact that $\mat{U}$ is unitary and $\mat{w}$ is diagonal: \begin{gather*} \mat{1} = \mat{U}\mat{U}^\dagger = \braket{\mat{P}|\mat{x}}\mat{w}\braket{\mat{x}|\mat{P}},\\ \delta_{ij} = \sum_{a}P_{i}^*(x_a)w_aP_{j}(x_a). \end{gather*} \begin{gather*} \mat{1} = \mat{U}^\dagger\mat{U} = \mat{w}^{1/2}\braket{\mat{x}|\mat{P}}\braket{\mat{P}|\mat{x}}\mat{w}^{1/2},\\ \delta_{ab} = w_{a}^{1/2}\sum_{n}P_{n}(x_a)P_{n}^*(x_b)w_{b}^{1/2} = w_{a}\sum_{n}P_{n}(x_a)P_{n}^*(x_b). \end{gather*} ::: ::::: :::{prf:theorem} DVR Bases -- braket notation Consider a pseudo-spectral bases of $N$ orthonormal functions $\ket{P_n}$, abscissa $\ket{x_a}$, and quadrature weights $w_a$. Let \begin{gather*} [\mat{G}]_{an} = \braket{x_a|P_n}, \qquad [\mat{w}]_{ab} = w_{a}\delta_{ab}, \qquad \mat{U} = \mat{G}^\dagger\mat{w}^{1/2}. \end{gather*} The following 5 statements are *equivalent*: 1. $\mat{G}\mat{G}^\dagger = \mat{w}^{-1}$ is diagonal. 2. $\mat{U}^\dagger = \mat{U}$ is unitary. 3. Discrete orthogonality: \begin{gather*} \sum_{a=0}^{N-1} w_a P_{i}^*(x_a)P_{j}(x_a) = \delta_{ij}, \qquad \sum_{a=0}^{N-1} w_a \braket{P_i|x_a}\braket{x_a|P_j} = \delta_{ij}, \qquad \Tr \ket{P_j}\bra{P_i} = \delta_{ij}, \qquad \end{gather*} 4. Discrete completeness: \begin{gather*} \sum_{n=0}^{N-1} w_a P_{n}^*(x_a)P_{n}(x_b) = \delta_{ab} \end{gather*} 5. Discrete $\delta$-property (locality): \begin{gather*} \phi_{n}(x_a) = \frac{\delta_{na}}{\sqrt{w_{n}}, \qquad \braket{\phi_m|\phi_n} = \delta_{mn}. \end{gather*} ::: The simplest path to a [DVR][] bases is from a set of $N$ orthonormal polynomial functions $P_{n}(x)$ and an $N$-point Gaussian quadrature rule -- see e.g. {cite}`Schneider:2004`. The more general approach, however, is to consider projected delta-functions {cite}`LCCMP:2002`, so we will start here. ### As Projected Delta Functions A general framework for [DVR][] bases {cite}`LCCMP:2002` starts with a projection operator $\op{P}$ which projects a Hilbert space onto states of low momentum. The [DVR][] basis is formed from projected delta-functions $\Delta_{n}(x) = \braket{x|\op{P}|x_n}$ with abscissa $x_m$ chosen so that they sit in the nodes: :::{margin} The last identify follows from the fact that $\op{P}^\dagger \op{P} = \op{P}^2 = \op{P}$. ::: \begin{gather*} \Delta_{n}(x_m) = K_{m}\delta_{mn} = \braket{\Delta_{n}|\Delta_{m}}. \end{gather*} From these, we define the orthonormal basis $\ket{F_m}$: \begin{gather*} \ket{F_{m}} = \frac{1}{\sqrt{K}_{n}}\ket{\Delta_{n}}, \qquad \braket{F_{m}|F_{n}} = \delta_{mn}, \qquad \braket{x|F_{n}} = F_{n}(x). \end{gather*} Note that this implies the following locality property \begin{gather*} \braket{x_m|F_n} = F_n(x_m) = \frac{\Delta_n(x_m)}{\sqrt{K_n}} = \sqrt{K_{n}}\delta_{mn}. \end{gather*} Thus, if such a basis and set of abscissa can be found, then for functions that live in the space spanned by the basis, the coefficients of the expansion can be found by simply evaluating the functions at the abscissa: \begin{gather*} \ket{f} = \sum_{n} \ket{F_n}\underbrace{f_n}_{\frac{f(x_n)}{\sqrt{K_{n}}}}, \qquad f(x_m) = \sum_{n} F_n(x_m)f_n = \sqrt{K_{m}}f_m. \end{gather*} :::{admonition} Example: Sinc-function Basis Consider for example, 1-D where we have \begin{gather*} \op{1} = \int_{-\infty}^{\infty}\!\!\!\!\d{x} \ket{x}\!\!\bra{x} = \int_{-\infty}^{\infty}\frac{\d{k}}{2\pi} \ket{k}\!\!\bra{k},\\ \braket{x|k} = e^{\I k x}, \qquad \braket{x|y} = \delta(x-y), \qquad \braket{k|q} = (2\pi)\delta(k-q). \end{gather*} Here we have the position basis vectors $\ket{x}$ (denoted by variables $x$, $y$, $z$, etc.) and momentum basis vectors $\ket{k}$ (denoted by variables $k$, $q$, $p$, etc.) In the position basis, $\ket{x}$ are Dirac delta-functions, while $\ket{k}$ are plane waves. We can now form the projection operator onto states with $\abs{k} \leq K$: \begin{gather*} \op{P} = \int_{-K}^{K}\frac{\d{k}}{2\pi} \ket{k}\!\!\bra{k}. \end{gather*} We use this to form the projected delta-functions: \begin{gather*} \Delta_{n}(x) = \braket{x|\Delta_n} \braket{x|\op{P}|x_n} = \int_{-K}^{K}\frac{\d{k}}{2\pi} \braket{x|k}\!\!\braket{k|x_n} = \int_{-K}^{K}\frac{\d{k}}{2\pi} e^{\I k(x_n - x)}\\ = \frac{e^{\I K(x_n - x)} - e^{-\I K(x_n - x)}} {2\pi\I(x_n - x)} = \frac{K}{\pi}\underbrace{\frac{\sin\Bigl(K(x_n - x)\Bigr)} {K(x_n - x)}}_{\sinc\Bigr(K(x-x_n)\Bigr)}. \end{gather*} Let's start with an abscissa $x_0 = 0$ so that \begin{gather*} \Delta_0(x) = \frac{K}{\pi}\sinc(Kx), \qquad x_{n} = \frac{\pi n}{K}, \end{gather*} where the abscissa are chosen to be the nodes. One can easily verify that locality property is satisfied, and use the fact that $\op{P}^\dagger \op{P} = \op{P}^2 = \op{P}$ to show that \begin{gather*} \braket{\Delta_{m}|\Delta_{n}} = \braket{x_m|\Delta_n} = \Delta_n(x_m) = K_{m} \delta_{mn}. \end{gather*} Note that this says that the projected delta-functions are orthogonal iff they are zero at all other abscissa. Evaluating the non-zero overlap, we find the normalization factors $K_{m} = \Delta_{m}(x_m) = K/\pi$, and can construct the orthonormal basis states $\ket{F_n}$: \begin{gather*} \ket{F_{n}} = \frac{1}{\sqrt{K_n}} \ket{\Delta_n} = \sqrt{\frac{\pi}{K}} \ket{\Delta_n}, \qquad \braket{x|F_{n}} = \sqrt{\frac{K}{\pi}}\sinc\Bigl(K(x - x_n)\bigr),\\ x_n = \frac{\pi n}{K}, \qquad \braket{F_m|F_n} = \delta_{mn}. \end{gather*} ::: :::{doit} Construct the [DVR][] basis for periodic function $x\in [0, L)$. Hint: Note that periodicity $f(x+L) = f(x)$ requires a specific form of $k$ for the basis states, so the projection operator becomes a sum rather than an integral. As a check, your results should be a linear combination of the sines and cosines the provide a basis for periodic functions (a.k.a. Fourier series). ::: To use this basis, one needs to compute the kinetic energy operator (Laplacian). This can be approximated using the locality property, e.g.: \begin{gather*} K_{mn} = -\braket{F_{m}|\op{\nabla}^2|F_{n}} \approx -F''_{n}(x_m). \end{gather*} A nice property is that, for many physically cases, the matrix for the potential can be well approximated by a diagonal matrix \begin{gather*} V_{mn} = \braket{F_m|V(\op{x})|F_n} \approx V(x_m)\delta_{mn}. \end{gather*} See {cite}`LCCMP:2002` and references therein for the conditions under which this approximation is good. :::{admonition} Example: Sinc-function Basis Cont. \begin{gather*} K_{mn} = -\braket{F_{m}|\op{\nabla}^2|F_{n}} = \braket{F'_{m}|F'_{n}}\\ = \frac{K}{\pi}\int_{-\infty}^{\infty}\d{x}\; \sinc'\Bigl(K(x-x_m)\Bigr)\sinc'\Bigl(K(x-x_n)\Bigr)\\ = \frac{K^2}{\pi}\int_{-\infty}^{\infty}\d{\theta}\; \sinc'(\theta)\sinc'\bigl(\theta - \pi(n-m)\bigr). \end{gather*} These integrals are a bit nasty, but if we are acting on functions well represented in the basis, such that $\op{\nabla}^2\ket{f}$ is also well represented in the basis, then we can use the locality property: \begin{gather*} K_{mn} \approx \frac{-F''_{n}(x_m)}{\sqrt{K_{m}}} \end{gather*} \begin{gather*} \sinc'(z) = \frac{z\cos{z} - \sin{z}}{z^2},\\ \sinc''(z) %= \frac{-z^2\sin{z} - 2z\cos{z} + 2\sin{z}}{z^3}\\ = \sinc{z}\left(\frac{2}{z^2}-1\right) - \frac{2\cos{z}}{z^2}. \end{gather*} If we consider $n \neq m$, then the first term will vanish \begin{gather*} K_{m\neq n} \approx \frac{2\cos\Bigl(K(x_m - x_n)\Bigr)}{(x_m - x_n)^2} = \frac{2K^2\cos\bigl(\pi(m-n)\bigr)}{\pi^2(m-n)^2} = \frac{2K^2(-1)^{n-m}}{\pi^2(m-n)^2}. \end{gather*} When $n=m$ we must consider all the terms and take the limit using l'Hopital's rule: \begin{gather*} \lim_{z\rightarrow 0} \sinc''(z) = \frac{-1}{3} K_{nn} = \frac{K^2}{3}. \end{gather*} Thus, \begin{gather*} K_{mn} \approx K^2 \begin{cases} \dfrac{1}{3} & n = m,\\ \dfrac{2(-1)^{n-m}}{\pi^2(m-n)^2}, & n \neq m. \end{cases} \end{gather*} It turns out that this is exactly the matrix you would get by computing the overlap integrals. ::: As a quick check, we consider the spectrum of a Harmonic oscillator with exact energies \begin{gather*} E_n = \hbar \omega (n + \tfrac{1}{2}). \end{gather*} If we want 20 states, then we must choose a basis where $\hbar^2 K^2/2m > \hbar\omega (20.5)$ giving an estimate of $K > 6.4/a_0$ where $a_0 = \sqrt{\hbar/m\omega}$ is the oscillator length. We also need to make sure that the basis is large enough in position space. This state will have turning point $R = \sqrt{2 E_{20}/m\omega^2} \approx 6.4a_0$ which means we need $n \pi / K > R$. A bit of playing shows that $Ka_0 = 10$ works well, which requires $n > 6.4 K/\pi \approx 20$. Again, playing, we find that $n = 30$ works well. Thus, we can get 20 states to machine precision with about 61 basis vectors. ```{code-cell} def get_err_HO(Ka=10.0, N=30): hbar = m = w = 1 a0 = np.sqrt(hbar/m/w) K = Ka/a0 n = np.arange(-N, N+1) x = np.pi * n / Ka V = np.diag(m*(w*x)**2/2) n_, m_ = n[None, :], n[:, None] D2 = K**2 * np.ma.divide(2*(-1.0)**(n_ - m_)/np.pi**2, (m_-n_)**2).filled(1/3) K = hbar**2/2/m * D2 H = K + V assert np.allclose(H, H.T) E = np.linalg.eigvalsh(H) exact = hbar*w * (np.arange(len(n)) + 0.5) err = abs(E/exact - 1) return err fig, ax = plt.subplots(figsize=(4,3)) ax.semilogy(get_err_HO(Ka=10.0, N=32)) ax.set(xlabel="$n$", ylabel="Energy error (rel)"); ``` ### From Orthogonal Polynomials If the basis functions $P_n(x)$ are orthogonal polynomials of degree $n$, then we have an alternative and simpler path to a [DVR][] basis. This section follows the ideas in {cite}`Schneider:2004` but using a consistent notation with the previous section. Start with a basis and $N$-point [Gaussian quadrature][] rule accurate for polynomials of degree $2N-1$: \begin{gather*} \braket{P_m|P_n} = \int w(x)\d{x}\; P_{n}^*(x)P_{m}(x) = \delta_{mn},\qquad \int w(x)\d{x}\; f(x) = \sum_{n} w_n f(x_n). \end{gather*} Specifically, the quadrature rule is valid for $f(x) = P_{n}^*(x)P_{m}(x)$ where $n+m < 2N$. The [DVR][] basis is then \begin{gather*} \ket{F_n} = \sum_{m}\ket{P_m}c_{mn}, \qquad F_n(x_m) = \sqrt{K_n}\delta_{mn}, \qquad \braket{F_{m}|F_{n}} = \delta_{mn}. \end{gather*} For this to hold, the orthonormality and quadrature rules imply \begin{gather*} c_{mn} = \braket{P_{m}|F_{n}} = \sum_{a} w_{a}P_{m}^*(x_a)F_{n}(x_a) %= \sum_{a} w_{a}P_{m}^*(x_a)\sqrt{K_{n}}\delta_{na} = w_{n}P_{m}^*(x_{n})\sqrt{K_{n}}\\ \braket{F_{m}|F_{n}} = \sum_{a} w_aF_m^*(x_a)F_n(x_a) = K_{n} w_{n}\delta_{mn}. \end{gather*} Hence \begin{gather*} K_{n} = \frac{1}{w_{n}}, \qquad c_{mn} = \sqrt{w_{n}}P_{m}^*(x_n),\qquad F_{n}(x) = \sqrt{w_{n}}\sum_{m}P_{m}^*(x_n)P_{m}(x). \end{gather*} Since $F_{n}(x)$ is a polynomial of degree $N-1$ with $N-1$ roots $\{x_{m\neq n}\}$, we can express it as a [Lagrange polynomial][]: \begin{gather*} F_{n}(x) = \frac{1}{\sqrt{w_n}}\prod_{m}\frac{x-x_{m}}{x_{n}-x_{m}}. \end{gather*} [Lagrange polynomial]: ### Some Technical Details Let's consider trying to construct a [DVR][] basis from a general set of $N$ orthogonal functions $\ket{P_n}$ and abscissa $x_n$: \begin{gather*} \ket{F_i} = \sum_{n}\ket{P_n}U_{ni}, \qquad \braket{x_m|F_i} = \sqrt{K_m}\delta_{mi}, \qquad \braket{F_a|F_b} = \delta_{ab}. \end{gather*} For these to be consistent we need \begin{gather*} \sum_{n}P_{n}(x_m)U_{ni} = \sqrt{K_m}\delta_{mn}. \end{gather*} In matrix notation, we need \begin{gather*} [\mat{M}]_{mn} = P_{n}(x_m), \qquad [\mat{U}]_{ni} = U_{ni},\\ \mat{M}\mat{U} = \mat{D} = \diag(\sqrt{K}), \qquad \mat{U}^\dagger\mat{U} = \mat{1}. \end{gather*} Thus, consistency requires that the matrix $\mat{M}$ can be diagonalized by right-multiplying by a unitary matrix $\mat{U}$. :::{margin} This is true in finite-dimensional spaces, but might need some care in infinite dimensional spaces since there $\mat{c}^\dagger\mat{c} = \mat{1}$ might not imply that $\mat{c}\mat{c}^\dagger = \mat{1}$. ::: This requires that \begin{gather*} \mat{M}\mat{U}\mat{U}^\dagger\mat{M}^\dagger = \mat{M}\mat{M}^\dagger = \abs{\mat{D}}^2 = \diag(K). \end{gather*} This means that the abscissa must be chosen so that the rows of $\mat{M}$. These are the vectors $P_n(x)$ consisting of the basis functions evaluated at each colocation point $x=x_{i}$. It is not clear we can always do this: We have $N$-degrees of freedom -- the locations of the abscissa $x_{n}$ -- but must satisfy $N(N-1)/2$ conditions for the upper triangular portion of $\mat{M}\mat{M}^{\dagger}$ to be zero. For orthonormal polynomials, this orthogonality follows from the [Christoffel–Darboux formula][] \begin{gather*} \sum_{n=0}^{N-1} \frac{P_{n}(x)P_{n}(y)}{\braket{P_n|P_n}} = \sum_{n=0}^{N-1} \braket{x|P_n}\braket{P_n|y} = A_{N}\frac{P_{N}(x)P_{N-1}(y) - P_{N-1}(x)P_{N}(y)}{x-y} \end{gather*} if we choose the abscissa to be the roots of $P_{N}(x_n)=0$. Other functions satisfy a similar identity, including the Bessel functions {cite}`Tygert:2006`, \begin{gather*} \sum_{k=1}^{\infty} 2(\nu+k)J_{\nu+k}(x)J_{\nu+k}(y) = \frac{xy}{x-y}\Bigl(J_{\nu+1}(x)J_{\nu}(y) - J_{\nu}(x)J_{\nu+1}(y)\Bigr)\\ 2(\nu+1)\sum_{k=1}^{\infty} 2(\nu+2k)J_{\nu+2k}(x)J_{\nu+2k}(y) = \frac{x^2y^2}{x^2-y^2}\Bigl(J_{\nu+2}(x)J_{\nu}(y) - J_{\nu}(x)J_{\nu+2}(y)\Bigr), \end{gather*} and Airy functions [Christoffel–Darboux formula]: ## Basis Functions We start with basis functions -- typically some set of [orthogonal polynomials][] $P_n(x)$, but these can be more general, including Fourier bases (sines and cosines). Most convenient are functions orthogonal with respect to some integration weight $w(x)$: \begin{gather*} \braket{P_{m}|P_{n}} = \int \d{x}\; w(x) P_{m}^*(x) P_{n}(x) = \delta_{mn}. \end{gather*} ::::::{tab-set} :::::{tab-item} Bessel The Bessel function basis is suitable for the radial coordinate in cylindrical and spherical geometries. close relationship with Fourier series gives rise to many interesting properties of this basis. \begin{gather*} P_{n}(x) = J_{\nu}(nx), \qquad P_{0}(x) = \frac{T_0(x)}{\sqrt{\pi}}, \qquad P_{n>0}(x) = \sqrt{\frac{2}{\pi}}T_{n}(x),\\ w(x) = \frac{1}{\sqrt{1-x^2}}, \qquad \braket{P_{m}|P_{n}} = \int_{-1}^{1} \d{x}\; \frac{P_{m}(x)P_{n}(x)}{\sqrt{1-x^2}} = \delta_{mn}. \end{gather*} ::::: :::::{tab-item} Fourier The Fourier basis consists of plane waves, orthogonal on a periodic interval: \begin{gather*} P_{n}(x) = e^{\I k_{n} x}, \qquad k_{n} = \frac{2\pi n}{L}, \qquad w(x) = \frac{1}{L}\\ \braket{k_{m}|k_{p}} = \frac{1}{L}\int_{0}^{L} \d{x}\; e^{-\I k_{m}x}e^{\I k_{n}x} = \delta_{mn}. \end{gather*} ::::: :::::{tab-item} Chebyshev The Chebyshev basis is a transformed cosine basis to the interval $x\in [-1, 1]$. The close relationship with Fourier series gives rise to many interesting properties of this basis. \begin{gather*} T_{n}(x) = \cos(n\cos^{-1} x), \qquad P_{0}(x) = \frac{T_0(x)}{\sqrt{\pi}}, \qquad P_{n>0}(x) = \sqrt{\frac{2}{\pi}}T_{n}(x),\\ w(x) = \frac{1}{\sqrt{1-x^2}}, \qquad \braket{P_{m}|P_{n}} = \int_{-1}^{1} \d{x}\; \frac{P_{m}(x)P_{n}(x)}{\sqrt{1-x^2}} = \delta_{mn}. \end{gather*} ::::: :::::: The idea is to express your functions in terms of this basis: \begin{gather*} f(x) = \sum_{n} a_n P_n(x). \end{gather*} Then the properties of the basis functions -- such as their derivatives -- can be used to calculate properties of $f(x)$. ::::::{admonition} Do It! Use a Fourier Basis to do the first example in {cite}`Boyd:1989`. In §1.2, {cite}`Boyd:1989` solves the following equation using a low-order polynomial expansion: \begin{gather*} u_{,xx} - (x^6 + 3x^2)u = 0, \qquad u(-1) = u(1) = 1. \end{gather*} Read the solution there, then repeat that analysis using a Fourier basis using everything you learn, including appropriately chosen collocation points, and considering the discussion about $a_1 = 0$. :::{admonition} {solution} As Boyd points out, the solutions must be even (prove this!), hence we need only consider a cosine series. The boundary conditions can be enforced with a series of the form \begin{gather*} u(x) = 1 + \sum_{n=0}^{\infty} a_n\cos(k_n x), \qquad k_n = \bigl(n-\tfrac{1}{2}\bigr)\pi. \end{gather*} Considering two terms $a_0$ and $a_1$, we choose the Gauss-Chebyshev collocation points: the roots of $\phi_{1}(x) = \cos(k_1 x) = \cos(\pi x/2) = 0$: \begin{gather*} k_0 = \frac{\pi}{2}, \qquad k_1 = \frac{3\pi}{2}, \\ x_{0} = \frac{1}{2}, \qquad x_{1} = 1. \end{gather*} Plugging the solution into the equation gives us the residual function, which we set to zero at each of the collocation points \begin{gather*} -(x_{i}^6 + 3x_{i}^2) + \sum_{n}\Bigl( -k_n^2 - (x_{i}^6 + 3x_{i}^2) \Bigr)a_n\cos(k_n x_{i}) = 0\\ -\left( \frac{\pi^2}{4} + \Bigl(\frac{1}{2^6} + \frac{3}{2^3}\Bigr) \right)a_0 -\left( \frac{\pi^2}{4} + \Bigl(\frac{1}{2^6} + \frac{3}{2^3}\Bigr) \right)(a_0\cos\tfrac{\pi = \Bigl(\frac{1}{2^6} + \frac{3}{2^3}\Bigr),\\ -\left( \frac{9\pi^2}{4} + \Bigl(1 + 3\Bigr) \right)a_1 = \Bigl(1+3\Bigr),\\ \begin{pmatrix} \frac{\pi^2}{4} + \frac{7}{8}\\ \end{pmatrix} \begin{pmatrix} a_0\\ a_1 \end{pmatrix} = \begin{pmatrix} \frac{7}{8}\\ 4 \end{pmatrix}N \end{gather*} ::: :::::: ```{code-cell} import sympy def get_spectral_approx(N, show=False): an = sympy.symbols(f'a:{N}', real=True) kn = tuple(sympy.pi * (2*n + 1)/2 for n in sympy.S(np.arange(N))) xi = tuple((2*n+1)/(2*N+1) for n in sympy.S(np.arange(N))) print(xi) #xi = tuple(sympy.cos(sympy.pi * (N-n)/2/N) for n in sympy.S(np.arange(N))) x = sympy.symbols('x') u = 1 + sum(a_ * sympy.cos(k_ * x) for a_, k_ in zip(an, kn)) Eqs = [(u.diff(x,x) - (x**6 + 3*x**2)*u).subs([(x, x_)]).evalf() for x_ in xi] sol = sympy.solve(Eqs, an) u_x = u.subs(sol) if show: display(sol) display(u.evalf()) u_ = sympy.lambdify(x, u_x, "numpy") return u_ fig, ax = plt.subplots() x = np.linspace(-1, 1) u_exact = np.exp((x**4-1)/4) ax.plot(x, u_exact, label="Exact: $u(x) = e^{(x^4-1)/4}$") for N in [2, 3, 4, 5, 6]: u = get_spectral_approx(N=N) ax.plot(x, u(x), label=f"Spectral: ${N=}$") fig.legend() ``` ## Numerical Example: Chebyshev Basis We will explore these more in {ref}`sec:ChebyshevBasis`, but for now, we give a couple of simplex examples using the Chebyshev polynomials $T_{n}(x)$, which form a basis for functions on the interval $x\in [-1, 1]$: \begin{gather*} T_{n}(x) = \cos(n\cos^{-1}x). \end{gather*} ### Example 1: Derivatives We start by using the basis the approximate a function $f(x)$ and then to compute its derivative. Let's consider $f(x) = e^{x}$: to determine the coefficients $a_{n}$ we define a **residual** which we minimize. There are many ways of doing this -- for example: \begin{gather*} R(\vect{a}) = \int_{-1}^{1}\frac{\d{x}}{\sqrt{1-x^2}} \Bigl(f(x) - \sum_{n}a_n T_n(x)\Bigr)^2. \end{gather*} To better work with this, we introduce the vectors $\ket{f}$ (the function) an $\ket{f_{\vect{a}}}$ (the approximation): \begin{gather*} \braket{f|g} = \int_{-1}^{1}\frac{\d{x}}{\sqrt{1-x^2}} f^*(x) g(x),\\ \ket{f_{\vect{a}}} = \sum_{n}\ket{T_n}a_{n} = \ket{\vect{T}}\cdot\vect{a},\\ \braket{x|f} = f(x), \qquad \braket{x|t_{\vect{x}}} = \sum_{n}a_{n}\braket{x|T_{n}} = \sum_{n}a_nT_{n}(x). \end{gather*} The residue, can then be expressed: \begin{gather*} R(\vect{a}) = \braket{f|f} - 2\braket{f|f_{\vect{a}}} + \braket{f_{\vect{a}}|f_{\vect{a}}}. \end{gather*} Approximate the ground state energy and wavefunction for a particle in a box with a Gaussian barrier: \begin{gather*} V(x) = \begin{cases} V_0 e^{-x^2/2\sigma^2} & -1 < x < 1,\\ \infty & \text{otherwise}. \end{cases} \end{gather*} ## DVR Bases :::{margin} Slightly more generally, we can define the inner-product in terms of a [Lebesque-Stieltjes integral][] with a measure $\alpha(x)$, which is any non-decreasing function on the real numbers $\alpha: \mathbb{R} \mapsto \mathbb{R}$: \begin{gather*} \braket{P_m|P_n} = \int P_m^*(x)P_n(x) \d{\alpha(x)}. \end{gather*} The relationship to the **weight function** $w(x)$ is $\d{\alpha(x)} = w(x)\d{x}$: i.e. $w(x) = \alpha'(x)$. For example, the [Lebesque-Stieltjes][Lebesque-Stieltjes integral] with $\alpha(x) = \Theta(x)$ rigorously defines what we mean by a Dirac delta function weight $w(x) = \delta(x)$. ::: To construct a [DVR][] basis, one needs in addition to the basis, a quadrature rule. Thus, we need: 1. A set of $N$ orthonormal basis functions $P_{n}(x)$. 2. An $N$-point [Gaussian quadrature][] rule exact for all linear and quadratic products of the basis functions. As before, orthonormality is defined with respect to an inner product \begin{gather*} \braket{P_m|P_n} = \int\d{x}\; w(x)P_m^*(x)P_n(x) = \delta_{mn}, \end{gather*} and the $N$-point quadrature rule consists of $N$ points $x_i$ and weights $w_i$ such that \begin{gather*} \int \d{x}\; w(x)f(x) = \sum_{i} w_i f(x_i). \end{gather*} For a DVR basis, we need this quadrature formula to be exact for: 1. The functions: $P_n(x)$. 2. Products: $P_m^*(x)P_n(x)$. 3. Products with an additional factor of $x$: $xP_m^*(x)P_n(x)$. *(This is not strictly needed, but ensures that keeping the potential diagonal works reasonably well.)* This will be the case for the [orthogonal polynomials][], which is why we use the notation $P_n(x)$, but also holds for the Fourier basis and some other sets with appropriately chosen quadrature rules. From these ingredients, we construct a set of **coordinate basis functions** $F_n(x)$ that are **local** in the sense that $F_n(x_m)$ vanishes unless $m=n$, and that are orthonormal with respect to the coordinate $x$: \begin{gather*} F_m(x) = \sum_{n} c_{mn}P_n(x), \qquad F_m(x_n) = f_m\delta_{mn}, \\ \braket{F_m|F_n} = \delta_{mn}, \qquad \braket{F_m|x|F_n} = x_m\delta_{mn}. \end{gather*} Using the orthonormality and evaluating the integrals with the quadrature rule, these properties give the relationship: \begin{gather*} c_{mn} = \int \d{x}\; w(x) P_{n}^*(x)F_{m}(x) = \sum_{i} w_i P_{n}^*(x_i)F_{m}(x_i) = w_m P_{n}^*(x_m)f_m,\\ \delta_{mn} = \int \d{x}\; w(x) F_m^*(x)F_n(x) = \sum_{i} w_i F_m^*(x_i)F_n(x_i) = w_m \abs{f_m}^2. \end{gather*} :::{admonition} Other Notations and the Cardinal Functions :class: dropdown To compare with the literature, note that {cite}`Schneider:2004` uses the notation $u_n(x) = \sqrt{w(x)}F_n(x)$ so that expectation values are simple \begin{gather*} \braket{u_m|A(x)|u_n} = \int \d{x}\; u_m^*(x) A(x) u_n(x). \end{gather*} This is convenient for numerical work. In the notation of {cite}`LCCMP:2002`, $K_m = 1/w_m$. They use the following three functions, which differ only by normalization: \begin{gather*} \Delta_n(x) = \sqrt{K_n}F_n(x) = K_n L_n(x),\\ \frac{\braket{\Delta_m|\Delta_n}}{K_m} = \braket{F_m|F_n} = K_m\braket{L_m|L_n} = \delta_{mn}\\ \frac{\Delta_{m}(x_n)}{K_{m}} = \frac{F_{m}(x_n)}{\sqrt{K_m}} = L_{m}(x_n) = \delta_{mn}. \end{gather*} The functions $L_{m}(x)$ here are denoted $C_{m}(x)$ in {cite}`Boyd:1989` and are called **cardinal functions**. These act like a [Dirac delta function][] in the basis and can often be obtained by projection. For example, the cardinal functions for the Fourier basis can be constructed from (for even $N$) \begin{gather*} C_{m}(x) = \braket{x|C_{m}} = \frac{1}{N}\sum_{n}\braket{x|P_n}\braket{P_n|x_{m}} = \frac{1}{N}\sum_{n}P_{n}(x)P_{n}^*(x_m)\\ = \frac{1}{N}\sum_{n=-N/2}^{N/2-1}e^{\I k_{n}(x-x_m)} = \frac{1}{N}\sum_{n=-N/2}^{N/2-1}e^{2\pi \I n(x-x_m)/L}. \end{gather*} This is geometric series: \begin{gather*} S = \sum_{n=a}^{b} r^a = \frac{r^{b+1} - r^{a}}{1-r}, \qquad r = e^{2\pi \I (x-x_m)/L}, \qquad -a = b+1 = \frac{N}{2}. \end{gather*} Thus, \begin{gather*} C_{m}(x) = \frac{1}{N}\frac{e^{\pi N \I (x-x_m)/L} - e^{-\pi N \I (x-x_m)/L}} {1-e^{2\pi\I (x-x_m)/L}} = e^{-\pi\I (x-x_m)/L} \frac{\sin\frac{\pi N (x-x_m)}{L}}{N\sin\frac{\pi(x-x_m)}{L}}. \end{gather*} The phase factor is somewhat inconsequential. *(It results from the fact that, for even $N$, there is a mismatch of positive and negative frequencies, so that the complex plane-wave $P_{-N/2}(x)$ is unpaired. This issue can be removed, e.g., by shifting $k_n \rightarrow k_n + \pi/L$, or by choosing $N$ to be odd.)* What remains are the cardinal functions for a periodic sine basis. If we take the limit $L \rightarrow \infty$ holding $\d{x} = L/N$ fixed, then we recover the Whitaker cardinal functions (see the [Whittaker-Shannon interpolation formula][]) for the real line: \begin{gather*} \lim_{N, L \rightarrow \infty} C_{m}(x) = \frac{\sin\Big(k_c (x-x_m)\Bigr)}{k_c(x-x_m)} = \sinc\Big(k_c (x-x_m)\Bigr), \qquad k_c = \frac{\pi N}{L} = \frac{\pi}{\d{x}}. \end{gather*} ::: These can be solved, taking $f_m$ to be real: \begin{gather*} f_m = \frac{1}{\sqrt{w_m}}, \qquad c_{mn} = \sqrt{w_m}P_{n}^*(x_m), \qquad F_{m}(x) = \sqrt{w_m}\sum_{n}P_{n}^*(x_m)P_n(x). \end{gather*} If the last quadrature condition is satisfied, then we have $\braket{F_m|x|F_n} = x_m\delta_{mn}$, which gets to the essence of these basis: that the potential can be expressed as a diagonal matrix: \begin{gather*} [\mat{V}]_{mn} = \braket{F_m|V(x)|F_n} = V(x_m)\delta_{mn}. \end{gather*} To compute the kinetic energy, we must include factors of $\sqrt{w(x)}$ in the basis functions: \begin{gather*} u_n(x) = \sqrt{w(x)}F_{n}(x). \end{gather*} Then, the kinetic energy can be evaluated using integration by parts and the quadrature rule. See below for details. ### Orthogonal Polynomials :::{margin} The integral $\int \d{\alpha(x)}$ is a [Lebesque-Stieltjes integral] with measure $\alpha(x)$, which should be a non-decreasing function on the reals. More commonly, we write $\d{\alpha(x)} = w(x)\d{x}$ where $w(x) = \alpha'(x)$ is called the **weight function**, but this is less general/more ambiguous. E.g. $\alpha(x) = \Theta(x)$ rigorously defines what we mean by a Dirac delta function $w(x) = \delta(x)$. ::: One way of constructing a DVR basis is from a set of $N$ polynomials $P_n(x)$ (i.e. of maximum degree $N-1$) that are orthonormal with respect to a measure $\alpha(x)$: \begin{gather*} \braket{P_m|P_n} = \int P_m(x)P_n(x) \d{\alpha(x)}. \end{gather*} An $N$-point [Gaussian quadrature][] rule is a set of $N$ points $x_i$ and weights $w_i$ such that \begin{gather*} \int P(x)\d{\alpha} = \sum w_i f(x_i) \end{gather*} is exact for all polynomials of degree $2N-1$ or less. To be used in constructing the DVR basis, the quadrature rule must be exact for polynomials of order $2(N-1) + 1 = 2N-1$, thus, the [Gaussian quadrature][] provides an acceptable set of $N$ abscissa and weights for constructing the DVR basis using the formulae above. ## Aliasing and Conservation Laws :::{margin} This discussion remains true in 3D. In 1D there are infinitely many conserved quantities since the system is integrable. I expect there are pseudo-spectral formulations with as many conserved quantities as degrees of freedom. ::: In the continuum limit with constant external potential, the GPE conserves particle number, energy, and momentum: \begin{align*} \DeclareMathOperator{\Im}{Im} N &= \int \abs{\psi}^2\d{x}, & E &= \int \left(\frac{\hbar^2\abs{\vect{\nabla}\psi}^2}{2m} + \frac{g\abs{\psi}^4}{2} \right)\d{x}, & \vect{ȷ} &= \Im \int (\psi^\dagger\hbar\vect{\nabla}\psi)\d{x}. \end{align*} To see this, use the GPE \begin{gather*} \I \hbar \dot{\psi} = \underbrace{\Biggl(\frac{-\hbar^2\nabla^2}{2m} + g\abs{\psi}^2\Biggr)} _{\op{H}_{\psi}}\psi. \end{gather*} Consider the conservation of momentum: \begin{gather*} \dot{\vect{ȷ}} = \Re \int \Biggl( \psi^\dagger\vect{\nabla}(\op{H}_{\psi}\psi) - (\op{H}_{\psi}\psi)^\dagger\vect{\nabla}\psi \Biggr) \d{x}. \end{gather*} The contributions from the kinetic energies poses will cancel since derivatives commute, but the non-linear potential requires an application of the product rule: \begin{gather*} \dot{\vect{ȷ}} = g\Re \int \Biggl( \psi^\dagger\vect{\nabla}(\abs{\psi}^2\psi) - \abs{\psi}^2 \psi^\dagger\vect{\nabla}\psi \Biggr) \d{x} = g\Re \int \Biggl( \psi^\dagger\psi\vect{\nabla}(\abs{\psi}^2) + \abs{\psi}^2\psi^\dagger\vect{\nabla}\psi - \abs{\psi}^2 \psi^\dagger\vect{\nabla}\psi \Biggr) \d{x},\\ = g\Re \int \Biggl(\abs{\psi}^2\vect{\nabla}(\abs{\psi}^2)\Biggr) \d{x},\ \end{gather*} ```{code-cell} Nx = 64 Lx = 10.0 dx = Lx/Nx k = 2*np.pi * np.fft.fftfreq(Nx, dx) k_max = abs(k).max() x = np.arange(Nx) * dx - Lx//2 Dk = 1j*k Dk[Nx//2] = 0 D2k = -k**2 def P(f, kc): return np.fft.ifft(np.where(abs(k) [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]: [Gaussian quadrature]: