12. Orthogonal Polynomials

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 Assignment 6: Orthonormal Bases, 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:

\[\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 [Abramowitz and Stegun, 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 [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*}\]

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)$.

Do It! Derive this condition.

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*}\]

Do It! Derive this solution from $(wp)' = wq$.

\[\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 [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*}\]

Theorem 2 (Riesz-Fischer – Theorem 7.2.1 in [Hassani, 2013])

The space $\mathcal{L}^2_{W}(a, b)$ is complete.

Theorem 3 (Theorem 7.2.2 in [Hassani, 2013])

All complete inner product spaces with countable bases are isomorphic to $\mathcal{L}^2_{W}(a, b)$.

Theorem 4 (Stone-Weierstrass approximation – Theorem 7.2.3 in [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.

Do It! Find a formula for the weights $w_i$.

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 [Press et al., 2007]):

\[\begin{gather*} w_{i} = \frac{\braket{f_{N-1}|f_{N-1}}}{f_{N-1}(x_i)f_{N}'(x_i)}. \end{gather*}\]

# 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}")

np.__version__='2.3.5', sp.__version__='1.16.3'
       Numpy, Formula above
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: False, False
N= 14: True , False
N= 15: False, False
N= 16: False, False
N= 17: False, False
N= 18: False, False
N= 19: False, False

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 scipy.integrate.quad() and 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

B_1 = -1/2
B_2 = 1/6
B_3 = 0
B_4 = -1/30
B_5 = 0
B_6 = 1/42
B_7 = 0
B_8 = -1/30
B_9 = 0
B_10 = 5/66
B_11 = 0
B_12 = -691/2730
B_13 = 0
B_14 = 7/6
B_15 = 0
B_16 = -3617/510
B_17 = 0
B_18 = 43867/798
B_19 = 0
B_20 = -174611/330

\[\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*}\]

Do It! Prove this.

\[\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)$.#

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.

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*}\]

12. Orthogonal Polynomials

Contents

12. Orthogonal Polynomials#

Approximation#

Quadrature#

Bernoulli Numbers and Polynomials#

Mittag-Lefler Theorem and \(\zeta(2k)\).#