--- jupytext: formats: ipynb,md:myst,py:light text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.8 kernelspec: display_name: Python 3 language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os 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 inline import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None ``` (sec:OrthogonalPolynomials)= # 12. Orthogonal Polynomials Generically, orthogonal polynomials can be scaled to provide an orthonormal basis for functions on some interval $[a, b]$ with respect to some weighting function $W(x)$: \begin{gather*} \braket{f|g} = \int_{a}^{b}f^*(x)g(x)\, W(x)\d{x}. \end{gather*} As demonstrated in {ref}`sec:Assignment6`, one can always construct these from the basis polynomials $\ket{x^{n}}$ using the [Gram-Schmidt process][]. Let the coefficients be called $k_n^{(i)}$: \begin{gather*} f_n(x) = k_n^{(n)}x^n + k_{n}^{(n-1)}x^{n-1}\cdots + k_{n}^{(0)}. \end{gather*} Some examples, with applications to quantum mechanics, are: :::{margin} Note: most classical orthogonal polynomials have a conventional scaling that makes them not orthonormal. For those listed, the conventional scalings are \begin{gather*} P_n(1) = P^*_{n}(1) = T_n(1) = 1, \\ L_n(x): \qquad k_n = \frac{(-1)^n}{n!},\\ H_n(x): \qquad k_n = 2^{n}. \end{gather*} See {cite}`Abramowitz:1965` for details. ::: \begin{align*} &\text{Legendre:} &&P_{n}(x), & x &\in [-1, 1], & W(x) &= 1, &&\text{(Spherical Harmonics)},\\ &\text{Shifted Legendre:} &&P^{*}_{n}(x), & x &\in [0, 1], & W(x) &= 1, &&\text{},\\ * [ ] &\text{Chebyshev (1st kind):} & &T_{n}(x), & x &\in [-1, 1], & W(x) &= \frac{1}{\sqrt{1-x^2}}, &&\text{(Fourier/Spectral Methods)},\\ &\text{Hermite:} &&H_{n}(x), & x &\in (-\infty, \infty), & W(x) &= e^{-x^2}, &&\text{(Harmonic Oscillator)},\\ &\text{Laguerre:} &&L_{n}(x), & x &\in [0, \infty), & W(x) &= e^{-x}, &&\text{(Hydrogen Atom)}. \end{align*} While programmatically straightforward, this process is tedious, and does not readily lead to easy-to-work-with expressions. Instead, one often prefers [Rodrigues formulas][] or recurrence formulas. From a practical perspective, most of the properties are well tabulated in resources like {cite}`Abramowitz:1965`. In particular, we have the following definitions: \begin{gather*} p(x) f_n'' + q(x)f_n' + \lambda_n f_n = 0, \tag{Sturm-Liouville ODE}\\ p(x) = \alpha x^2 + \beta x + \gamma, \qquad q(x) = \mu x + \nu, \qquad \lambda_n = -n(n-1)\alpha -n \nu,\\ f_n = \frac{1}{e_n W(x)}\diff[n]{}{x}\Bigl(W(x) [p(x)]^{n}\Bigr), \tag{Rodrigues formula} \end{gather*} Rodrigues formula is valid if $p(x)$ is at most quadratic, and $q(x)$ is linear. The weighting function comes from the self-adjointness condition \begin{gather*} W(x) = \frac{1}{p(x)}\exp\left(\int^{x}\frac{q(x)}{p(x)}\d{x}\right). \end{gather*} Additional, we have the following recurrence formula (see {cite}`Hassani:2013` Eq. (7.12)): \begin{gather*} f_{n+1}(x) = (a_n x + b_n)f_{n}(x) - c_{n}f_{n-1}(x), \tag{Recurrence Relation}\\ f_{n}(x) = k^{(n)}_{n}x^n + k^{(n-1)}_{n}x^{n-1} + \dots,\\ a_n = \frac{k^{(n+1)}_{n+1}}{k^{(n)}_{n}}, \qquad b_n = a_n\left( \frac{k^{(n)}_{n+1}}{k^{(n+1)}_{n+1}} - \frac{k^{(n-1)}_{n}}{k^{(n)}_{n}} \right),\qquad c_n = \frac{\braket{f_n|f_n}}{\braket{f_{n-1}|f_{n-1}}}\frac{a_n}{a_{n-1}}. \end{gather*} :::::{admonition} Example: Hermite Polynomials The (physicist's) [Hermite polynomials][] $H_n(x)$ satisfy: \begin{gather*} H_n(x) = (-1)^{n}e^{x^2}\diff[n]{}{x}e^{-x^2}, \tag{Rodrigues formula}\\ H_n''(x) - 2x H_n'(x) + 2n H_n(x) = 0, \tag{Sturm-Liouville}\\ H_{n+1}(x) = 2xH_n(x) - 2nH_{n-1}(x), \qquad H_0(x) = 1, \qquad H_1(x) = 2x, \tag{Recurrence} \end{gather*} Alternatively, $\psi_n(x) = e^{-x^2/2}H_n(x)$ satisfies the more familiar harmonic oscillator equation: \begin{gather*} -\frac{\psi''_n(x)}{2} + \frac{x^2}{2}\psi_n(x) = (n+\tfrac{1}{2})\psi_n(x). \end{gather*} Depending on your needs, one or more of these definitions may be useful. ::::: You should recall that one of the key properties of self-adjoint (hermitian) operators is that they admit a complete orthonormal basis. Thus it is natural to expect that orthogonal polynomials might arise as the eigenfunctions for some sort of hermitian ODE such as a [Sturm-Liouville][] problem. This is the basis for [Rodrigues formulas][]. Recall that a general second-order ODE of the form \begin{gather*} \mathcal{L} y = p(x)y'' + q(x)y' + \lambda y = 0 \end{gather*} is self-adjoint if $p'(x) = q(x)$. :::{admonition} Do It! Derive this condition. :class: dropdown Recall that, to be self-adjoint, we must have \begin{gather*} \braket{f|\mathcal{L}|g} = \braket{g|\mathcal{L}|f}\\ \int f \mathcal{L}g\d{x} = \int g \mathcal{L} f\d{x}. \end{gather*} Starting with \begin{gather*} \braket{f|\mathcal{L}|g} \end{gather*} and integrating by parts, we obtain (up to boundary terms, which we neglect for now) the condition $fpg'' + fqg' = gpf'' + gqf'$ \begin{align*} fpg'' + fqg' &= (pf)''g - (fq)'g = (p'f + pf')'g - (f'q + fq')g\\ &= (p''f + 2p'f' + pf'')g - (f'q + fq')g = gpf'' + gqf',\\ 0 &= (p'' - q')fg + 2(p' - q)f'g = (p' - q)'fg + 2(p' - q)f'g. \end{align*} This is clearly satisfied if $p' = q$ when we can write \begin{gather*} (py')' + \lambda y = 0. \end{gather*} ::: We can thus make a general equation self-adjoint by multiplying by a weight $W$ chosen such that \begin{gather*} (Wp)' = Wq, \qquad W(x) = \frac{1}{p(x)}\exp\left(\int^x \frac{q(x)}{p(x)}\d{x}\right). \end{gather*} :::{admonition} Do It! Derive this solution from $(wp)' = wq$. :class: dropdown \begin{gather*} (wp)' = w'p + wp' = wq, \qquad w' = w\frac{q-p'}{p}\\ \int\frac{\d{w}}{w} = \ln w = \int^{x}\frac{q-p'}{p}\d{x} = \int^x\frac{q}{p}\d{x} - \int^p\frac{\d{p}}{p} = \int\frac{q}{p}\d{x} - \ln p. \end{gather*} ::: With this weight, we have \begin{gather*} (wpy')' + \lambda w y = 0 \end{gather*} and the eigenfunctions can be chosen to be orthonormal with the inner product \begin{gather*} \braket{f|g}_{w} = \int f^*(x) g(x)\, W(x)\d{x}. \end{gather*} It is not clear that these eigenfunctions are polynomials, and in general, there will be solutions to the differential equation that are not -- these are the generalized functions associated with the orthogonal polynomials. ## Approximation Orthogonal polynomials are extremely useful for approximations. Part of this follows from two theorems discussed in {cite}`Hassani:2013`. These refer to the space of square-integrable functions $f(x)$: \begin{gather*} \mathcal{L}^2_{W}(a, b) = \left\{f(x) \text{ where } \int_{a}^{b} \abs{f(x)}^2\, W(x)\d{x} = \braket{f|f} \text{ is finite}\right\}. \end{gather*} :::{prf:theorem} Riesz-Fischer -- Theorem 7.2.1 in {cite}`Hassani:2013` :label: thm:RieszFisher The space $\mathcal{L}^2_{W}(a, b)$ is complete. ::: :::{prf:theorem} Theorem 7.2.2 in {cite}`Hassani:2013` All complete inner product spaces with countable bases are isomorphic to $\mathcal{L}^2_{W}(a, b)$. ::: :::{prf:theorem} Stone-Weierstrass approximation -- Theorem 7.2.3 in {cite}`Hassani:2013` The set of monomials $\{x^k\}_{k=0}^{\infty}$ forms a basis for $\mathcal{L}^2_{W}(a, b)$. ::: This means that we can approximate any function in $\mathcal{L}^2_{W}(a, b)$ with polynomials: \begin{gather*} f(x) \approx \sum_{k=0}^{N} a_k x^k. \end{gather*} To do this with brute force, we might evaluate the function at a set of $N$ abscissa $[\vect{f}]_{n} = f(x_n)$, then try to interplate: \begin{gather*} f(x_n) = \sum_{k=0}^{N} a_k x_{n}^k = \sum_{k}V_{nk}a_{k}\\ \vect{f} = \mat{V}\vect{a}. \end{gather*} The matrix $\mat{V}$ is called a [Vandermonde matrix][] and in can be very poorly conditioned making it hard to solve for $\vect{a}$ numerically. In contrast, once you have an orthonormal basis, once can find the coefficients trivially (and stably) by simply computing the overlap integrals. ## Quadrature Another useful property of orthogonal polynomials is [numerical integration][] or **quadrature**. In general, given $N$ abscissa $\{x_i\}$ one can define $N$ weights $w_i$ such that \begin{gather*} \int_{a}^{b} f(x)\, W(x)\d{x} = \sum_{i=0}^{N-1} w_i f(x_i) \end{gather*} is exact for polynomials $f(x)$ of degree $N-1$ or less. :::{admonition} Do It! Find a formula for the weights $w_i$. :class: dropdown Solve the following linear system of $N$ equations for $n \in \{0, 1, \dots, N-1\}: \begin{gather*} \int_{a}^{b} x^n\, W(x)\d{x} = \sum_{i=0}^{N-1} w_i x_i^{n}. \end{gather*} This gives a set of weights $w_i$ that are exact for the monomials $x^{n}$. Since integration is a linear operation, these are therefore exact for all linear combinations -- hence for all polynomials. Note that solving for $w_i$ requires inverting the following [Vandermonde matrix][] \begin{gather*} [\mat{V}]_{ni} = x_{i}^{n}. \end{gather*} For equally spaced abscissa, this matrix becomes ill-conditioned, so this is not a good numerical strategy. Orthogonal polynomials provide a much better approach. ::: If we can also choose the abscissa, then we can use the additional freedom to make the quadrature rule exact for polynomials up to order $2N-1$. Such a set of abscissa $\{x_i\}$ and weights $\{w_i\}$ is called a [Gaussian quadrature][] rule. Once you have a set of orthogonal polynomials for a give weight function $W(x)$ and interval $[a, b]$, the $N$-point [Gaussian quadrature] rule can be obtained by using the $N$ roots of $f_n(x_i) = 0$ with weights (see {cite}`PTVF:2007`): \begin{gather*} w_{i} = \frac{\braket{f_{N-1}|f_{N-1}}}{f_{N-1}(x_i)f_{N}'(x_i)}. \end{gather*} ```{code-cell} # Example: Hermite Polynomials import scipy as sp from scipy.integrate import quad print(f"{np.__version__=}, {sp.__version__=}") def W(x): return np.exp(-x**2) def braket(f, g): return quad(lambda x: f(x)*g(x)*W(x), -np.inf, np.inf)[0] print(" Numpy, Formula above") for N in range(1,20): # Monic polynomials f_N = np.polynomial.hermite.Hermite([0]*(N) + [1/(2)**(N)]) f_N1 = np.polynomial.hermite.Hermite([0]*(N-1) + [1/(2)**(N-1)]) x = f_N.roots() w = braket(f_N1, f_N1) / f_N1(x) / f_N.deriv()(x) x_, w_ = np.polynomial.hermite.hermgauss(N) assert np.allclose(x, x_) assert np.allclose(w, w_) M = 2*N ai = np.array([quad(lambda x: x**n * W(x), -np.inf, np.inf)[0] for n in range(M)]) aq = np.array([sum(sorted(x**n * w)) for n in range(M)]) aq_ = np.array([sum(sorted(x_**n * w_)) for n in range(M)]) print(f"{N=: 3}: {np.allclose(ai, aq_)!s:5}, {np.allclose(ai, aq)!s:5}") ``` Note, however, that one must be very careful with round-off error. The results above differ, even on the same hardware but between different versions of NumPy and SciPy. *(I have not traced through where the difference is coming from, but all results are different, so there is something affecting both {func}`scipy.integrate.quad` and {func}`numpy.polynomial.hermite.hermgauss`.)* ``` np.__version__='1.26.4', sp.__version__='1.15.2' Numpy, Formula N= 1: True , True N= 2: True , True N= 3: True , True N= 4: True , True N= 5: True , True N= 6: True , True N= 7: True , True N= 8: True , True N= 9: True , True N= 10: True , True N= 11: True , False N= 12: True , False N= 13: True , False N= 14: True , False N= 15: False, False N= 16: True , False N= 17: True , False N= 18: False, False N= 19: False, False ``` # Bernoulli Numbers and Polynomials The [Bernoulli polynomials][] $B_n(x)$ defined by the generating function \begin{gather*} g(x, t) = \frac{t e^{xt}}{e^t - 1} = \sum_{n=0}^{\infty} B_n(x) \frac{t^n}{n!},\\ B_{n}(x) = \sum_{k=0}^{n} \binom{n}{k}B_{n-k}x^{k}, \end{gather*} where $B_{n} = B_n(0)$ are the [Bernoulli numbers][] ```{code-cell} :tags: [hide-input, margin] from fractions import Fraction from scipy.integrate import quad from scipy.special import factorial def f(theta, n, r=5.0): z = r * np.exp(1j*theta) dz_dtheta = 1j*z return factorial(n) * z/(np.exp(z) - 1)/z**(n+1) * dz_dtheta / (2j*np.pi) def get_B(n, r=5.0): res = quad(f, 0, 2*np.pi, args=(n, r), complex_func=True)[0] assert np.allclose(res.imag, 0) return Fraction(res.real).limit_denominator(10000) #for n in [0, 1] + [2*n for n in range(1, 11)]: for n in range(1, 21): print(f"B_{n} = {get_B(n)}") ``` \begin{gather*} g(0,t) = \frac{t}{e^{t} - 1} = \sum_{n=0}^{\infty} \frac{B_n}{n!}t^{n},\\ B_n = \left.\diff[n]{}{t}\frac{t}{e^t-1}\right|_{t=0} = n! \oint \frac{z}{e^{z}-1}\frac{1}{z^{n+1}}\frac{\d{z}}{2\pi\I}. \end{gather*} In the last formula, the contour should be taken to be a small circle about $z=0$ within the radius of convergence for the series. Note that \begin{gather*} \partial_x g(x, t) = tg(x, t)\\ \sum_{n=0}^{\infty} B_n'(x)\frac{t^{n}}{n!} = \sum_{n=0}^{\infty} B_n(x)\frac{t^{n+1}}{n!} = \sum_{n=1}^{\infty} B_{n-1}(x)\frac{t^{n}}{(n-1)!} = \sum_{n=1}^{\infty} nB_{n-1}(x)\frac{t^{n}}{n!}, \end{gather*} which leads to the useful property \begin{gather*} B_0'(x) = 0, \qquad B_n'(x) = nB_{n-1}(x). \end{gather*} Expanding a few terms: \begin{gather*} g(0,t) = \frac{t}{1 + t + \frac{t^2}{2!} + \cdots - 1} = \frac{1}{1 + \frac{t}{2!} + \cdots} = 1 - \frac{t}{2!} + \cdots. \end{gather*} Thus \begin{gather*} B_0 = 1, \qquad B_1 = -\frac{1}{2}. \end{gather*} Subtracting off the linear term we find an even function: \begin{gather*} g(0,t) + \frac{t}{2} = \frac{t}{2}\frac{e^{t} + 1}{e^{t} - 1} = \frac{t}{2}\frac{e^{t/2} + e^{-t/2}}{e^{t/2} - e^{-t/2}} = \frac{t/2}{\tanh(t/2)}. \end{gather*} Thus, all the remaining terms are even: \begin{gather*} g(0,t) = \underbrace{- \frac{t}{2}}_{\frac{B_{1}}{1!}t} + \sum_{n \text{ even}} \frac{B_{n}}{n!}t^n, \qquad B_{3} = B_{5} = B_{7} = \dots = 0. \end{gather*} Some algebra gives \begin{gather*} B_{2} = \frac{1}{6}, \qquad B_{4} = \frac{-1}{30}, \qquad B_{6} = \frac{1}{42}, \qquad B_{8} = \frac{-1}{30}, \qquad \cdots. \end{gather*} These are useful in the following series: \begin{gather*} \cot(x) = \frac{2\I}{y}\sum_{n\text{ even}}\frac{B_{n}}{n!}y^{n} = \frac{1}{x}\sum_{n\text{ even}}\frac{(-1)^{n/2}2^{n}B_{n}}{n!}x^{n}, \qquad y = 2\I x. \end{gather*} The series for $\tan x$ follows from \begin{gather*} \tan x = \cot x - 2 \cot 2x. \end{gather*} :::{admonition} Do It! Prove this. :class: dropdown \begin{align*} \cot(x) = \frac{\cos(x)}{\sin(x)} &= \I\frac{e^{\I x} + e^{-\I x}}{e^{\I x} - e^{-\I x}}\\ &= \I\left(\frac{e^{y/2} + e^{-y/2}}{e^{y/2} - e^{-y/2}} + 1 - 1\right)\\ &= \I\left(1 + \frac{2e^{-y/2}}{e^{y/2} - e^{-y/2}}\right)\\ &= \I\left(1 + \frac{2}{e^{y} - 1}\right)\\ &= \frac{2\I}{y}\left(\frac{y}{2} + \frac{y}{e^{y} - 1}\right)\\ &= \frac{2\I}{y}\sum_{n\text{ even}}\frac{B_{n}y^{n}}{n!}. \end{align*} ::: ## Mittag-Lefler Theorem and $\zeta(2k)$. :::{margin} A [meromorphic function][] is a function that is complex analytic (differentiable) except at a set of isolated points $z_n$ which are poles. I.e., it is a function that is analytic everywhere except at some poles. Every [meromorphic function][] can be expressed as a ratio of [holomorphic function][]s -- those which are analytic (complex differentiable) everywhere. *(Both these terms refer to some domain, which in this case is the entire plane.)* ::: The Mittag-Lefler theorem states that, if a [meromorphic function][] $f(z)$ falls off fast enough $\lim_{\abs{z}\rightarrow \infty} f(z)/\abs{z} \rightarrow 0$, then \begin{gather*} f(z) = f(0) + \sum_{n=1}^{N} r_n\left(\frac{1}{z-z_n} + \frac{1}{z_n}\right) \end{gather*} where $r_n$ are the residues of the respective poles. :::{admonition} Do It! Prove this. Hint: consider \begin{gather*} I_N = \oint_{C_N}\frac{f(w)}{w(w-z)}\frac{\d{w}}{2\pi \I} \end{gather*} for a contour $C_{N}$ enclosing the first $N$ poles (with smallest magnitude) of $f(z)$. ::: Noting that $\cot x$ is periodic with poles at $x_n = n \pi$ and residues (evaluated with [L'Hôpital's rule][]) \begin{gather*} r_n = \lim_{x\rightarrow x_n} (x-x_n)\cot(x) = \lim_{x\rightarrow x_n}\cos(x)\frac{x-x_n}{\sin(x)} = \cos{x_n}\lim_{x\rightarrow x_n}\frac{1}{\cos(x)} = 1. \end{gather*} \begin{gather*} \cot{z} = \frac{1}{z} + \sum_{n\neq 0}\left(\frac{1}{z-n\pi} + \frac{1}{n\pi}\right)\\ = \frac{1}{z} + \sum_{n>0}\left(\frac{1}{z-n\pi} + \frac{1}{z+n\pi}\right)\\ = \frac{1}{z} + \sum_{n>0}\frac{2z}{z^2-n^2\pi^2}\\ z\cot{z} = 1 + \sum_{n>0}\frac{2z^2}{z^2-n^2\pi^2}\\ \pi z\cot(\pi z) = 1 + \sum_{n>0}\frac{2\pi^2z^2}{\pi^2 z^2-n^2\pi^2} = 1 + \sum_{n>0}\frac{2z^2}{z^2-n^2} = 1 - \sum_{n>0}\frac{2z^2/n^2}{1-\frac{z^2}{n^2}}. \end{gather*} Expanding each term as a geometric series in $z^2$, we have \begin{gather*} \frac{1}{1-\frac{z^2}{n^2}} = \sum_{k=0}^{\infty} \left(\frac{z^2}{n^2}\right)^{k}\\ \pi z\cot(\pi z) = 1 - \sum_{n>0}\sum_{k=0}^{\infty}2\frac{z^2}{n^2}\frac{z^{2k}}{n^{2k}} = 1 - 2\sum_{k=1}^{\infty}z^{2k} \underbrace{\sum_{n=1}^{\infty}\frac{1}{n^{2k}}}_{\zeta(2k)}. \end{gather*} Hence \begin{gather*} 1 - 2\sum_{k=1}^{\infty}z^{2k}\zeta(2k) = \sum_{k=0}^{\infty}\frac{(-1)^{k}2^{2k}B_{2k}}{(2k)!}(\pi z)^{2k}. \end{gather*} Equating terms \begin{gather*} 1 = B_{0}, \qquad \zeta(2k) = \frac{(-1)^{k+1}B_{2k}(2\pi)^{2k}}{2(2k)!}. \end{gather*} [L'Hôpital's rule]: [meromorphic function]: [holomorphic function]: [Bernoulli polynomials]: [Bernoulli numbers]: [Gaussian quadrature]: [Vandermonde matrix]: [Gram-Schmidt process]: [Sturm-Liouville]: [Rodrigues formulas]: [Hermite polynomials]: [numerical integration]: