The Dirac Function: Series Approximation

狄拉克函数：级数逼近 苏剑林 · 2017-01-11 · Mathematics series approximation

The Weierstrass Theorem

Understanding the Dirac function as a limit of functions leads to a wealth of results, and moreover these results are not far removed from rigorous proofs. For example, define

$$\delta_n(x)=\left\{\begin{aligned}&\frac{(1-x^2)^n}{I_n},x\in[-1,1]\\ &0,\text{other cases}\end{aligned}\right.$$

where $I_n = \int_{-1}^1 (1-x^2)^n dx$, and it is not hard to show that

$$\delta(x)=\lim_{n\to\infty}\delta_n(x)$$

Then, for a continuous function $f(x)$ on $[a,b]$, we obtain

$$f(x)=\int_{-1}^1 f(y)\delta(x-y)dy = \lim_{n\to\infty}\int_{-1}^1 f(y)\delta_n(x-y) dy$$

Here $-1 < a < b < 1$, and we have "unrigorously" swapped the order of the integral and the limit — but that's not particularly important. What matters is the result: we can see that

$$P_n(x)=\int_{-1}^1 f(y)\delta_n(x-y) dy$$

is a degree-$2n$ polynomial in $x$, so the above equation shows that $f(x)$ is the limit of a sequence of degree-$2n$ polynomials! This leads to the famous "Weierstrass theorem":

Any continuous function on a closed interval can be uniformly approximated by polynomials.

Readers who insist on rigor will shake their heads: what kind of proof is this? Indeed, it isn't a rigorous proof — but how far is it from one? Readers interested in seeing a rigorous proof might look at Professor Qi Minyou's Revisiting Calculus, where they will find that the rigorous proof is really just a more careful discussion and bounding of the formulas above, and this discussion is not difficult. So although the above is entirely non-rigorous as a guide, it isn't far from a rigorous proof. All of this owes to our loose interpretation of $\delta(x)$ as the limit of $\delta_n(x)$, rather than treating it as a rigorous generalized function (distribution).

The process above actually reveals a general framework that lets us analyze what kind of function basis can be used to approximate a given function. Constructing a sequence of functions $\{\delta_n(x)\}$ that serves as a limit for $\delta(x)$ is not hard: first normalize it, then require that both tails of $x=0$ go to 0 as $n\to\infty$ — that's all there is to it. For example, here's another instance:

$$\delta_n(x)=\left\{\begin{aligned}&\frac{\cos^n x}{I_n},x\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right]\\ &0,\text{other cases}\end{aligned}\right.$$

where $I_n = \int_{-\pi/2}^{\pi/2} \cos^n x dx$, and again we have

$$\delta(x)=\lim_{n\to\infty}\delta_n(x)$$

What results can this limit give us? Following the same line of reasoning as with polynomials, suppose $- \frac{\pi}{2} < a < b < \frac{\pi}{2}$; then for a continuous function $f(x)$ on $[a,b]$ we find

$$\int_{-\pi/2}^{\pi/2} f(y)\delta_n(x-y) dy = \frac{1}{I_n} \int_{-\pi/2}^{\pi/2} f(y)\cos^n (x-y) dy $$

Expanding $\cos^n(x-y)$, we see that the above expression is actually a linear combination of $\sin kx, \cos kx, k=0,1,\dots,n$, so we also obtain:

Any continuous function on the closed interval $[-\pi/2, \pi/2]$ can be uniformly approximated by the sine-cosine series $\{\cos kx, \sin kx\}$.

This is also called the Weierstrass theorem. Likewise, the discussion above is not far from a rigorous proof.

Does something feel off here? $\cos x,\sin x$ has period $2\pi$, yet all we've approximated is the interval $[-\pi/2,\pi/2]$ of length $\pi$ — isn't that a bit wasteful? In fact, by slightly modifying the process above, one can prove:

For any $\epsilon > 0$, a continuous function on the closed interval $[-\pi+\epsilon, \pi-\epsilon]$ can be uniformly approximated by the sine-cosine series $\{\cos kx, \sin kx\}$.

However, this fails to hold if we replace $[-\pi, \pi]$ with — convergence can only be in measure or in the $L^2$-norm; as is well known, the Fourier basis does not give uniform convergence.

Orthonormal Bases

"Orthogonality" helps simplify problems: for geometric problems we like orthogonal coordinate systems, and for function problems we like orthogonal bases.

The two Weierstrass theorems above tell us that a continuous function on a closed interval can be uniformly approximated either by polynomials or by a sine-cosine series, and the process above actually also gives us a way to compute the approximation coefficients for each term. But from a practical computation standpoint, this approach is impractical, because going from an order-$n$ approximation to an order-$n+1$ approximation requires recomputing all the coefficients from scratch (by comparison, consider the Taylor series approximation: going from order $n$ to order $n+1$ only requires computing one additional term, $f^{(n+1)}(x)$ — very economical, but the Taylor series has much stronger requirements). This is quite wasteful, so we need to find a more efficient computational scheme.

We restrict our discussion to a closed interval $[a,b]$. For two functions $f(x),g(x)$, we define their inner product as

$$\langle f, g\rangle=\int_a^b f(x)g(x)dx$$

Suppose a sequence of functions $\{e_k (x)\},\, k=0,1,2,\dots$ forms a function basis satisfying

$$\langle e_m, e_n\rangle=\delta_{mn}$$

Then it is called an orthonormal basis. With this, we can describe the "optimal approximation" of a function.

Suppose we have a function $f(x)$ on the closed interval $[a,b]$, and we know that $f(x)$ can be uniformly approximated by the orthonormal basis $\{e_k (x)\},\, k=0,1,2,\dots$. Consider the finite partial sum

$$\sum_{k=0}^n \alpha_k e_k(x)$$

We want to use this finite series to optimally approximate $f(x)$, meaning we want to minimize the error, defined as

$$E=\int_a^b \left(f(x)-\sum_{k=0}^n \alpha_k e_k(x)\right)^2 dx$$

By "optimal" we mean adjusting $\alpha_0,\alpha_1,\dots,\alpha_n$ so that $E$ attains its minimum value. Since this is a finite series, we can directly take partial derivatives with respect to each $\alpha_k$ and solve, finding that when

$$\alpha_k = \int_a^b f(x)e_k (x)dx$$

$E$ is minimized. This gives us the optimal approximation. Notice that in this case, going from an order-$n$ approximation to an order-$n+1$ approximation only requires computing one additional term, $\alpha_{n+1} = \int_a^b f(x)e_{n+1} (x)dx$ — there's no need to recompute all the coefficients.

Note that if $\{e_k (x)\}$ is an arbitrary sequence of functions (not a basis, and not necessarily orthogonal), the same optimization procedure can still be carried out, but there's no guarantee that the resulting approximation actually converges to the original function. That is, it's possible that no matter how large you make $n$, the approximation error never decreases. However, if it has already been proven that $\{e_k (x)\}$ can uniformly approximate the original function, then we can guarantee that the error from the above optimization procedure tends to 0 — because these are already optimal solutions, and if the optimal solution's error doesn't tend to 0, how could it possibly converge uniformly to the original function?

Fourier Series

The question now is: how do we obtain an orthonormal basis? By the Weierstrass theorems, we already know that both the power functions $1,x,x^2,\dots$ and the sine-cosine functions $1,\sin x, \cos x, \sin 2x, \cos 2x, \dots$ can serve as bases; let's focus on the latter for now. On the interval $[-\pi, \pi]$, one can check that $1,\sin x, \cos x, \sin 2x, \cos 2x, \dots$ is orthogonal, i.e., for any natural number $m,n$,

$$\left\{\begin{aligned}&\int_{-\pi}^{\pi} \cos mx \cos nx = 0\,( m\neq n)\\ &\int_{-\pi}^{\pi} \sin mx \sin nx = 0\,( m\neq n)\\ &\int_{-\pi}^{\pi} \sin mx \cos nx = 0 \end{aligned}\right.$$

and also

$$\int_{-\pi}^{\pi} \sin^2 mx dx = \int_{-\pi}^{\pi} \cos^2 mx dx = \pi$$

Therefore

$$\sqrt{\frac{1}{\pi}},\sqrt{\frac{1}{\pi}}\sin x, \sqrt{\frac{1}{\pi}}\cos x, \sqrt{\frac{1}{\pi}}\sin 2x, \sqrt{\frac{1}{\pi}}\cos 2x, \dots$$

forms an orthonormal basis. This gives us

$$f(x) = \gamma + \sum_{n=1}^{\infty} \alpha_n \sin nx + \sum_{n=1}^{\infty} \beta_n \cos nx$$

where

$$\begin{aligned}\gamma=&\frac{1}{\pi}\int_{-\pi}^{\pi} f(x)dx,\\ \alpha_n=&\frac{1}{\pi}\int_{-\pi}^{\pi} f(x)\sin nx dx,\\ \beta_n=&\frac{1}{\pi}\int_{-\pi}^{\pi} f(x)\cos nx dx \end{aligned}$$

This is exactly the Fourier series. But since we're working on the interval $[-\pi,\pi]$, the Fourier series can only converge in measure or in the $L^2$-norm, not uniformly. And on the interval $[-\pi/2, \pi/2]$, $1,\sin x, \cos x, \sin 2x, \cos 2x, \dots$ is not orthogonal, so an orthogonalization step is needed.

Orthogonal Polynomials

We've already discussed the approximation problem for sine-cosine series, arriving at the Fourier series. So what about the other basis from the Weierstrass theorem — the power-function basis? How should we handle that? The natural orthogonality of the sine and cosine functions greatly simplifies the problem, but power functions lack this orthogonality, so we have no choice but to perform an orthogonalization procedure.

The best-known orthogonalization procedure is Gram-Schmidt orthogonalization, which is covered in any basic linear algebra course, so we won't repeat it here. Applying Gram-Schmidt orthogonalization to the basis $\{x^n\},\,n=0,1,2,\dots$ on the interval $[-1,1]$ yields:

$$\begin{array}{c|c} \hline n & P_n(x) \\ \hline 0 & 1 \\ 1 & x \\ 2 & \frac{1}{2}(3x^2-1)\\ 3 & \frac{1}{2}(5x^3-3x)\\ 4 & \frac{1}{8}(35x^4-30x^2+3)\\ \vdots & \vdots\\ \hline\end{array}$$

These polynomials $P_n(x)$ are called the "Legendre polynomials." They are orthogonal, but not yet normalized. After normalizing them we get $\hat{P}_n(x)$, and we can then write down a result analogous to the Fourier series:

$$f(x)=\sum_{n=0}^{\infty} \alpha_n \hat{P}_n (x)$$

where

$$\alpha_n = \int_{-1}^1 f(x)\hat{P}_n (x) dx$$

Other Series Approximations

Approximation by sine and cosine functions can also be viewed as approximation using the imaginary exponential $e^{ikx},\,k=0,\pm 1,\pm 2,\dots$ as a basis. A natural question is: can we instead approximate using the real exponential $e^{kx}$?

Clearly, based on our earlier experience, this comes down to whether we can construct a Dirac delta sequence out of exponential functions. Consider the hyperbolic cosine $\cosh x$: it's an even function that, near $x=0$, behaves like an upward-opening parabola. If we consider the interval $[-1, 1]$, we can define

$$\delta_n (x) = \frac{(\cosh 1 - \cosh x)^n}{I_n},\quad I_n = \int_{-1}^1 (\cosh 1 - \cosh x)^n dx$$

It's not hard to show that $\lim_{n\to\infty} \delta_n (x) = \delta (x)$, so we can consider the sequence of functions

$$\frac{1}{I_n}\int_{-1}^1 f(y)[\cosh 1 - \cosh (x-y)]^n dy $$

and use this sequence to approximate $f(x)$. Expanding it directly shows that it is a linear combination of $\cosh kx,\sinh kx, k = 0,1,2,\dots,n$, i.e., a linear combination of $e^{kx},\,k=0,\pm 1,\pm 2,\dots, \pm n$. So indeed, real exponentials can be used for approximation.

This raises an interesting question: since both imaginary and real exponentials can be used for approximation, why has only the imaginary-exponential version (the Fourier series) received widespread study? My guess is that it's mainly because the real-exponential version isn't as elegant or as practical. All these approximations operate on a finite interval; the Fourier series, however, corresponds to periodic functions, so studying it on one finite interval is equivalent to studying it everywhere. Real exponentials, by contrast, only apply within a single interval and offer no help in extending to the whole line — not to mention that they simultaneously involve $e^x$ and $e^{-x}$, which causes divergence at both $x\to+\infty$ and $x\to -\infty$, giving them no advantage at all.