Approximation theory is the field of analysis that explores how to approximate things. It does this in the broadest sense - it deals with approximating functions, constants, distances, and providing algorithms to approximate anything that could be approximated. It is likely one of the most useful pieces of pure mathematics, applications ranging from approximating Einstein’s equations for relativity to constructing computer graphics using the Bézier curves formed by the Bernstein polynomials.

This project is divided is divided into three main chapters: one on functional analysis, one on the approximation of $\pi$, and one on polynomial interpolation and its application to sequences.

The approximation of constants has existed since almost the dawn of mathematics, beginning to be necessary with the Greek and Indian discovery of the irrational following the Pythagorean theorem. However, the other fields of approximation theory are much more modern, the reason for this being that a precise notion and definition of a function is required for them, which was only first developed with the work of Bernoulli in the late 18th century.

Approximation theory then developed as a supposed solution to several problems, including Euler’s work on minimizing error on the maps of Russia in attempting converge the non Euclidian geometry of the Earth into a 2 dimensional rectangular map. Aproximation theory greatly simplified this process. A similar problem was investigated by Nash in the 1950s.

From this followed Chebyshev and Fourier’s work on the approximation of functions, culminating in the work of Bernstein in his proof of the Weierstrass approximation theorem, proving that any function can be approximated by a polynomial. Polynomial interpolation was first developed through the work of Newton in the late 17th century, and reached its culmination with the work of Lagrange in 1795.

Because approximation theory is so broad, it shows how mathematics has evolved, beginning with Archimedes’ geometric approximations, showing the development of calculus and Newton’s use of the ’divided difference’ as the derivative was developed, going through to Taylor series and integration.

Now if only it was still 2016, then we would have a square date today, since $$=5004225042016 = 5004^2)

Functional approximation

Introduction

Given a continuous function $f(x)$ defined over the closed interval $[a,b]$, how can we represent $f$ as another function $g$ such that the difference between the two functions is as small as possible? Can we construct a summation, or infinite series in terms of simple functions such that the resulting approximation is as close as possible to the target function? The theorem which proves that we can do this is known as the Weierstrass approximation theorem which states that given function $f(x)$ defined over closed interval $[a, b]$ and $\epsilon >0$, it is possible to construct polynomial $P_n(x)$ such that $|f(x) - P_n(x)| < \epsilon$, and therefore that any function can be approximated as closely as desired by a polynomial. The Fourier and Chebyshev series on which I will elaborate below use the method of infinite summation of easier functions to approximate more complicated functions, and have astoundingly fast convergence rate. Although the mathematics asserts that the perfect approximation is formed when the series is evaluated at infinity, if we wish to evaluate this in practice, we can truncate the series at a value $N$, forming a compromise between computational energy and accuracy.

Fourier analysis

Fourier analysis is the approximation of an often complicated function using only the sums of the sine and cosine functions, proposed by French mathematician Joseph Fourier, who in his classic 1822 paper "The analytical study of heat" showed that decomposing more complex functions into simpler ones made it much easier to solve differential equations. Specifically, Fourier was trying to solve the heat equation, which had been one of the biggest problems of the time. However, by simplifying the functions in the differential equation to just the sinusoidal functions, whose differential properties were well understood, Fourier managed to solve the heat equation. This method has become known as ’solving by Fourier series’. To demonstrate how the Fourier series works, consider a basic sine wave:

Now, consider the graph of $\sin(2x)$:

Now, consider this graph:

Note how dissimilar this graph is, on first examination, to the original two graphs. Given just this graph (and with no prior experience of Fourier analysis) you would assume that the graph is that of a strange and irregular function. However, this graph is just the graph of the function: $$f(x) = \sin(x) + \sin(2x)$$ This is the base for Fourier analysis. Given any function, how can we simplify it into just the sum of the sine and cosine functions?

This gives rise to two fields of approximation:

1. The Fourier series: Any continuous, differentiable oscillatory function can be (using the Fourier series) described as the sum of trigonometric or exponential functions. In order to ensure convergence, the Dirichlet conditions much be met - these being that the function must be periodic, absolutely integrable over its period, have a finite number of discontinuities and a finite number of maxima and minima.

2. The Fourier transform: Any continuous, differentiable function satisfying all the Dirichlet conditions except periodicity, can be described as a continuous integral of a complex exponential function.

Fourier Series

Let us first examine the Fourier series. What we want to achieve is an infinite series approximation for $f(x)$, in terms of the sine and cosine functions, weighted by index $n$ which has the effect of reducing the wavelength of the function and increasing accuracy and convergence of the approximation. We need to work out formulas for coefficients $a_0$, $a_n$ and $b_n$: $$f(x) = a_0 + \sum^{\infty}_{n=1}a_n \cos nx + \sum^{\infty}_{n=1}b_n \sin nx$$ In practice, when using the Fourier series to approximate a function, we will not evaluate the series to infinity, we will instead cut it off at a certain value of iterations, in a compromise between computational efficiency and convergence. Note that we do not need to have formula for $b_n$ $\sin 0 = 0$.

Deriving the coefficients

How can we find $a_n$ and $b_n$? Recall first from linear algebra that one of the most basic operations which can be applied to two vectors is the dot product. We refer to the dot product as a dot in the middle of the multiplication (eg. $\vec{a} \cdot \vec{b}$). Given two vectors $\vec{a}, \vec{b} \in \mathbb{R}^n$, their dot product is given by either $$\vec{a} \cdot \vec{b} = |a||b|\cos\theta,$$ where $|x|$ is the length of vector $x$ and $\theta$ is the angle between the two vectors or: $$\vec{a} \cdot \vec{b} = a_1b_1 + ... a_nb_n$$ The first formula has quite an elegant geometric interpretation - consider, for the sake of example, two vectors $\vec{a}$ and $\vec{b}$ in 2 dimensional space.

Consider the line formed by extending the line formed by connecting the tip of $\vec{b}$ and the origin. The dot product is formed by projecting $\vec{a}$ onto this line, and then multiplying the length of the projected $\vec{a}$ by the length of $\vec{b}$, the result is the dot product of the two vectors. This is why the cosine teresting property of this is that this property does not change whether $\vec{b}$ is being projected unto the line formed by $\vec{a}$ or vise versa, both cases give the dot product. One more property of this dot product is that if the two vectors are orthogonal ($\theta = 90$) the dot product evaluates to 0, whereas if they point in opposite directions, their dot product is negative. Let us furthermore consider how the second formula which states that if given two nth dimensional vectors, we can establish the dot product by simply summating the product of each element in both vectors (this will be of great use to us later) corresponds to the principle of projecting one vector unto another.

How can we apply this, to test whether two given functions $f$ and $g$ are orthogonal? We must represent our function as a vector, which we could do by forming a single row vector of infinitesimal increments, and then forming the dot product, which in the vector results in multiplying all the values in each row and summating this. This is essentially equivalent to the Riemann sum over $a$ and $b$, therefore the dot product in the vector of our function is equal to $$\int^b_a f(x)g(x)dx$$

Where $a$ and $b$ are the ends of the function.

If we apply this to the trigonometric functions we are working with, we can take their dot product such that $$\langle \cos nx, \cos kx \rangle = \int^\pi_{-\pi} (\cos nx)(\cos kx)\,dx$$ To evaluate this integral, we use the product-to-sum identity of the cosine function, which states that $$\cos mx \cos nx = \frac{1}{2}[\cos(m+n)x + \cos(m-n)x]$$ Our integral therefore simplifies to: $$\langle \cos nx, \cos kx \rangle = \frac{1}{2}\left[\int^\pi_{-\pi}\cos(m+n)x\,dx + \int^\pi_{-\pi}\cos(m-n)x\,dx\right]$$ Therefore, by the first Fundamental Theorem of Calculus, we can say that $$\langle \cos nx, \cos kx \rangle = \frac{1}{2}[\frac{1}{m+n}\sin(m+n)x \bigg\rvert^\pi_{-\pi} + \frac{1}{m-n}\sin(m-n)x \bigg\rvert^\pi_{-\pi}]$$ As the sine function crosses the $x$ axis for multiples of $\pi$, $\frac{1}{m+n}\sin(m+n)x \bigg\rvert^\pi_{-\pi} + \frac{1}{m-n}\sin(m-n)x \bigg\rvert^\pi_{-\pi}$ will evaluate to 0, therefore, for $m \neq k$, . $$\langle \cos nx, \cos kx \rangle = 0$$ For $m = k$, both terms become: $$\langle \cos nx, \cos kx \rangle = \frac{1}{2}\left[\int^\pi_{-\pi}\cos(2nx)\,dx + \int^\pi_{-\pi}\cos(0)\,dx\right]$$ The integral of $\cos(0)$ over any interval of $x$ is simply $2\pi$, while the integral of $\cos(2nx)$ over a full period of $2\pi$ evaluates to 0. Thus, when $m = k$, $$\langle \cos nx, \cos kx \rangle = \frac{1}{2}(0 + 2\pi) = \pi$$ Therefore: $$\langle \cos nx, \cos kx \rangle = \begin{cases} 0 & \text{if } m \neq k \\ \pi & \text{if } m = k \end{cases}$$ This means that multiples of cosine $m$ and $k$ are orthogonal in the dot product for $m \neq k$ and are $2\pi$ in the dot product for $m = k$. We consider whether this applies to the sine function by similarly taking $$\langle \cos nx, \cos kx \rangle \int^{\pi}_{-\pi} (\sin nx)(\sin kx)\,dx$$ we use the product-to-sum identities to simplify the integrand. The product-to-sum identities state that: $$\sin(\alpha)\sin(\beta) = \frac{1}{2}[\cos(\alpha - \beta) - \cos(\alpha + \beta)].$$ Using this identity, we rewrite the integrand as: $$\begin{aligned} \int_{-\pi}^{\pi} (\sin nx)(\sin kx)\,dx &= \int_{-\pi}^{\pi} \frac{1}{2}[\cos((n-k)x) - \cos((n+k)x)]\,dx \\ &= \frac{1}{2}\left[\int_{-\pi}^{\pi} \cos((n-k)x)\,dx - \int_{-\pi}^{\pi} \cos((n+k)x)\,dx\right]. \end{aligned}$$ Now, we integrate term by term. Since the integral of the cosine function over a full period is zero, we have, by the First Fundamental Theorem of Calculus: $$\int_{-\pi}^{\pi} \cos((n-k)x)\,dx = 0,$$ and similarly, $$\int_{-\pi}^{\pi} \cos((n+k)x)\,dx = 0.$$ We can therefore state that, for all $m \neq n$: $$\langle \sin nx \sin mx \rangle = 0$$ When $m = k$, the integral becomes: $$\langle \sin nx, \sin kx \rangle = \frac{1}{2}\left[\int_{-\pi}^{\pi} \cos(0)\,dx - \int_{-\pi}^{\pi} \cos(2nx)\,dx\right].$$ The integral of $\cos(0)$ over any interval of $x$ is $2\pi$, while the integral of $\cos(2nx)$ over a full period of $2\pi$ evaluates to 0. Thus, when $m = k$, $$\langle \sin nx, \sin kx \rangle = \frac{1}{2}(2\pi - 0) = \pi.$$ Therefore, similarly to cosine, in the dot product, sine is orthogonal for all multiples $m \neq k$ and returns $\pi$ for $m = k$ , such that: $$\langle \sin nx, \sin kx \rangle = \begin{cases} 0 & \text{if } m \neq k \\ \pi & \text{if } m = k \end{cases}$$ We finally consider whether the dot product of $\sin x$ and $\cos x$ (weighted by $n$ and $k$ respectively) demonstrates an orthogonality between the sinusoidal functions over their period. We consider the integral: $$\langle \sin nx, \cos kx \rangle = \int^{\pi}_{-\pi} (\sin nx)(\cos kx)\, dx$$ Using the product-to-sum identities, we rewrite the integrand as: $$\sin(nx) \cos(kx) = \frac{1}{2} \left[ \sin((n+k)x) + \sin((n-k)x) \right]$$

Thus, the integral becomes: $$\langle \sin nx, \cos kx \rangle = \frac{1}{2} \left[ \int_{-\pi}^{\pi} \sin((n+k)x) \, dx + \int_{-\pi}^{\pi} \sin((n-k)x) \, dx \right]$$

Evaluating each integral separately, by the first Fundamental Theorem of Calculus, we find: $$\int_{-\pi}^{\pi} \sin((n+k)x) \, dx = \int_{-\pi}^{\pi} \sin((n-k)x) \, dx = 0$$

This holds true for all $n$ and $k$, including the case where $n = k$.

Since both integrals evaluate to zero, we conclude: $$\int_{-\pi}^{\pi} \sin(nx) \cos(kx) \, dx = 0$$ And therefore that, for all $n$ and $k$, the functions $\sin x$ and $\cos x$ are orthogonal by the dot product

This is another reason why the sinusoidal functions were so useful to Fourier as he was solving the heat equation, as well as the fact that the differential properties of the individual functions are well understood, when they are multiplied we can also very easily simplify differential equations thus formed.

We can use this orthogonality principle to attempt to separate out $a_n$ and $b_n$ from our original approximation for $f(x)$: $$\int ^\pi _{-\pi} f(x) \cos (kx)dx = \int ^\pi _{-\pi} (\sum^{\infty}_{n=0}a_n \cos nx + \sum^{\infty}_{n=1}b_n \sin nx)\cos (kx)dx$$ The integrals demonstrated earlier show that $\int^\pi_{-\pi} f(x) \cos(kx)\, dx = 0$ for all $n \neq k$. When $n = k$, we have: $$\int ^\pi _{-\pi} f(x) \cos (kx)dx = \int ^\pi _{-\pi} a_k (\cos kx)^2\,dx = a_k\pi$$ Therefore we have discovered a formula for $a_k$, namely that: $$a_k = \frac{1}{\pi}\int^\pi_{-\pi} f(x)\cos(kx)\,dx$$ for $k \neq 0$. The key to finding $a_0$ is to realise that, to ensure that the Fourier series is accurate, we need to include a constant term that represents the average value of the function over one period. The average value of a function is given by $$\frac{1}{2L}\int^L_{-L} f(x) \,dx$$

By integrating the function over its entire period the formula captures the total area under the curve of the function over one complete cycle. Dividing the integral by the length of the period scales the integral result to represent an average over the entire period. This normalization ensures that the average value is independent of the period’s length, focusing solely on the function’s behavior within one complete cycle. The result of this computation gives the average height of the function over one period. In graphical terms, if you were to draw the function’s graph over one period, the height of the rectangle whose area equals the integral would represent the average value of the function over that period. Therefore, we can say that:

$$a_0 = \frac{1}{2\pi}\int^\pi_{-\pi}f(x)\,dx$$ To achieve coefficient $b_n$, we multiply by $\sin kx$ and integrate, such that $$\int ^\pi _{-\pi} f(x) \sin (kx)\,dx = \int ^\pi _{-\pi} (\sum^{\infty}_{n=0}a_n \cos nx + \sum^{\infty}_{n=1}b_n \sin nx)\sin (kx)\,dx$$ Similarly, using the orthogonality principle, the integral will cancel to 0, for all $n \neq k$. When $n = k$, we have: $$\int ^\pi _{-\pi} f(x) \sin (kx)\,dx = \int ^\pi _{-\pi} b_k (\sin kx)^2\,dx = b_k\pi$$ which yields $$b_n = \frac{1}{2\pi}\int^\pi_{-\pi}f(x)\sin(kx)dx.$$ We therefore have the complete Fourier series: $$f(x) = \frac{1}{4\pi}\int^\pi_{-\pi}f(x)dx+\sum^{\infty}_{n=1}\frac{1}{\pi}\int^\pi_{-\pi} f(x)\cos(nx)dx \cos nx + \sum^{\infty}_{n=1}\frac{1}{2\pi}\int^\pi_{-\pi}f(x)\sin(nx)dx \sin nx$$ Note that, by beginning index n a t 1 for both series, we add a separate formula for $a_0$. Note also that this must be divided by two so that $a_0$ be distributed equally between sine and cosine.

An important thing to note whilst approximating a function using the Fourier series is which of the functions sine and cosine should be used. Given an even function (where $f(x) = f(-x)$) we use just the even cosine function, whereas given an odd function (where $f(-x) = -f(x)$) we use just the odd sine function. In all other cases we use both.

Diagram reconstructed for this web edition; original image not supplied.

Dirac Delta Function

Let us now use an example to illustrate the use of the Fourier series, namely the Dirac delta function. The Dirac delta function, denoted by $\delta(x)$, is shown in Figure 1.1 above. The Dirac delta function is very interesting, because it is equal to zero at all values $x$ except at $x=0$. This function is very good for describing oscillatory behaviour - for example in an article I recently wrote for the most recent Axiom edition, I described how the periodic synchronization of fireflies can be modelled by showing that the probability $P_N(T_b)$ of all $N$ fireflies flashing at time time $T_b$ must converge to the Dirac delta function.

We can loosely define this function as

$$\delta (x) \simeq \begin{cases} +\infty, & \text{if } x = 0 \\ 0, & \text{if } x \neq 0 \end{cases}$$

However, this piecewise definition is not very elegant. It is better to define the delta function in terms of its properties - namely that $$\int_{-\infty}^{\infty} \delta(x) \, dx = 1$$ Note that this function is even, $\delta(x) = \delta(-x)$ therefore when using the Fourier series to approximate it we will not use sines, and only cosines. Therefore, if we approximate the delta function between $-\pi$ and $\pi$, we have: $$\delta(x) = \sum^\infty_{n=0} a_n \cos nx$$ We also have the previously derived formula for $a_0$ $$a_0 = \frac{1}{2 \pi}\int^\pi_{-\pi}\delta(x)dx$$ We also know from the earlier given property that the integral of the delta function is always equal to one. Therefore: $$a_0 = \frac{1}{2\pi}$$ Now let us also evaluate $a_k$. $$a_k = \frac{1}{\pi}\int^\pi_{-\pi} \delta(x)\cos(kx)dx$$ We can again use the previously named property. Furthermore, we must remember that the delta function returns 0 at all values except x = 0, therefore it isn’t really necessary for us to evaluate the integral as the multiplication will go to zero at all points except x = 0. In other words - in the multiplication, x = 0 is the only point which is interesting enough to consider. This point is very convenient as cos(0) = 1. Therefore, $$a_k = \frac{1}{\pi}$$ This formula for the constant tells us something amazing and very specific to the Dirac delta function - the constants remain constant regardless of $k$. All the frequencies are the same. Normally, we would see a decay in values a $a_k$, however, the singular delta function has constant coefficients throughout. The Fourier series for the Dirac delta function is therefore: $$\delta(x) = \frac{1}{2\pi} + \frac{1}{\pi}(\sum^\infty_{n=1}\cos nx)$$

Complex Fourier exponential series

It is important for the sake of the Fourier transform which follows, that our series can be rewritten in terms of only a complex coefficient $c_n$. We can do this because, if we have the set of complex exponentials $e^{inx}$ where $n \in \mathbb{Z}$, we can say that these represent the base of the sine and cosine functions, where $n$ is the frequency of the wave. This is encapsulated by Euler’s formula that $$e^{ix} = \cos x + i\sin x.$$ Euler’s formula is a beautiful mathematical identity that relates the exponential function, sine function, and cosine function in a single equation:

$$e^{i\theta} = \cos(\theta) + i\sin(\theta)$$

This formula is quite elegant because it combines these fundamental mathematical functions in a concise and profound way. This formula is derived using power series expansions. It can be shown that the exponential function $e^x$ can be represented as a power series:

$$e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots$$

Furthermore, the sine function can also be represented as a power series:

$$\sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots$$

Similarly, the cosine function can be represented as a power series:

$$\cos(x) = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + \cdots$$

Now, let’s consider the expression $e^{ix}$:

$$e^{ix} = 1 + ix + \frac{(ix)^2}{2!} + \frac{(ix)^3}{3!} + \frac{(ix)^4}{4!} + \cdots$$

Simplify each term:

$$e^{ix} = 1 + ix - \frac{x^2}{2!} - i\frac{x^3}{3!} + \frac{x^4}{4!} + i\frac{x^5}{5!} - \cdots$$

Now, let’s group the real and imaginary parts:

$$e^{ix} = (1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \cdots) + i(x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots)$$

Notice now that the real part is exactly the power series expansion of $\cos(x)$, and the imaginary part is exactly the power series expansion of $\sin(x)$. Thus, we have:

$$e^{ix} = \cos(x) + i\sin(x)$$ This formula can be made even more beautiful if we let $x = \pi$, yielding $$e^{i\pi} + 1 = 0$$ which beautifully combines all the fundamental mathematical constants. Euler’s formula neatly represents the link between the exponential and sinusoidal functions, showing that we can rewrite our approimation in terms of just a complex coefficient $c_n$ and the exponential function. Therefore, we can approximate $f(x)$ as $$f(x) = \sum^{\infty}_{n=-\infty} c_n e^{inx},$$ Furthermore, $e^{inx}$ shares a similar orthogonality principle with the sinusoidal functions. This can be shown by evaluating $$\langle e^{inx}e^{-ikx} = \int^\pi_{-\pi}e^{inx}e^{-ikx}$$ To solvce this, we can use Euler’s formula, derived earlier by the power series. We express $e^{inx}$ and $e^{-ikx}$ by Euler’s formula:

$$e^{inx} = \cos(nx) + i\sin(nx)$$

$$e^{-ikx} = \cos(kx) - i\sin(kx)$$

Substituting these expressions into the integral, we get:

$$\int_{-\pi}^{\pi} (\cos(nx) + i\sin(nx))(\cos(kx) - i\sin(kx)) \, dx$$

Expanding this expression, we get:

$$\begin{align*} &\int_{-\pi}^{\pi} \cos(nx)\cos(kx) \, dx - i\int_{-\pi}^{\pi} \sin(nx)\cos(kx) \, dx \\ &\quad + i\int_{-\pi}^{\pi} \cos(nx)\sin(kx) \, dx + \int_{-\pi}^{\pi} \sin(nx)\sin(kx) \, dx \end{align*}$$

Now, since $\cos$ is an even function and $\sin$ is an odd function, the integrals involving products of $\cos$ and $\sin$ over a symmetric interval like $[-\pi, \pi]$ will all evaluate to zero. So we are left with:

$$\int_{-\pi}^{\pi} \cos(nx)\cos(kx) \, dx + \int_{-\pi}^{\pi} \sin(nx)\sin(kx) \, dx$$

Using the orthogonality property of sine and cosine functions, which states that:

$$\int_{-\pi}^{\pi} \cos(mx)\cos(nx) \, dx = \begin{cases} \pi & \text{if } m = n \\ 0 & \text{if } m \neq n \end{cases}$$

and similarly for $\sin$, we can evaluate the remaining integrals:

$$\int_{-\pi}^{\pi} \cos(nx)\cos(kx) \, dx = \begin{cases} \pi & \text{if } n = k \\ 0 & \text{if } n \neq k \end{cases}$$

$$\int_{-\pi}^{\pi} \sin(nx)\sin(kx) \, dx = \begin{cases} \pi & \text{if } n = k \\ 0 & \text{if } n \neq k \end{cases}$$

So, the integral becomes:

$$\pi(\delta_{nk} + \delta_{nk}) = \begin{cases} 2\pi & \text{if } n = k \\ 0 & \text{if } n \neq k \end{cases}$$

where $\delta_{nk}$ is the Kronecker delta, defined as:

$$\delta_{nk} = \begin{cases} 1 & \text{if } n = k \\ 0 & \text{if } n \neq k \end{cases}$$

Therefore,

$$\int_{-\pi}^{\pi} e^{inx}e^{-ikx} \, dx = \begin{cases} 2\pi & \text{if } n = k \\ 0 & \text{if } n \neq k \end{cases}$$ Demonstrating that, similar to the sinusoidal functions, the exponential functions are orthogonal by the dot product with different weights $n$ and $k$ for $n \neq k$, whereas for $n = k$, the dot product returns 0. Using this, we can derive the complex coefficient $c_n$, by multiplying the equation by $e^{-ikx}$ and integrating over the closed interval $[-\pi, \pi]$ to yield $$\int^\pi_{-\pi} f(x) e^{-ikx}\, dx = \int^\pi_{-\pi} (\sum^\infty_{-\infty} c_n e^{inx})e^{-ikx}\,dx$$ $$= \sum^\infty_{-\infty} c_n \int^\pi_{-\pi} e^{inx}e^{-ikx}\, dx$$ Therefore, for all $n \neq k$ the integral is 0. When $n = k$, it simplifies to $2\pi$, therefore we can write that $$\int^\pi_{-\pi} f(x) e^{-ikx}\, dx = 2\pi c_k$$ From this, we can say that $$c_k = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(x) e^{-ikx} \, dx$$ We can therefore write the complex exponential form of the Fourier series as $$f(x) = \sum^\infty_{-\infty}\frac{1}{2\pi} \int_{-\pi}^{\pi} f(x) e^{-inx} \, dx \,e^{inx}$$

Fourier Transform

We will use the complex form of the Fourier series, which I wrote to be $$f(x) = \sum^\infty_{-\infty}\frac{1}{2\pi} \int_{-\pi}^{\pi} f(x) e^{-inx} \, dx \,e^{inx}$$ If we extrapolate from the sinusoidal functions with base period $\pi$ into an arbitrary (must satisfy all the Dirichlet conditions except that of periodicity) function $f(x)$ of period $L$, the whole formula can be rewritten in terms of the period L as: $$f(x) = \sum^{\infty}_{n=-\infty} c_n e^{\frac{inx\pi}{L}}$$ If we let $$c_n = \frac{1}{2L}\int^L_{-L} f(x) e^{\frac{-in\pi x}{L}}dx$$ where L is the interval on which we set the function. Now, in order to be able to take the Fourier transform, we want to know what happens when we take the limit $$\lim_{L \to \infty} \sum^{\infty}_{n=-\infty} c_n e^{\frac{inx\pi}{L}}.$$ This means we no longer will be considering a fixed period, but the entire real numbers. We can try to make the equations slightly less convoluted by defining $$k_n = \frac{\pi n}{L}$$ It is important to note that when we summate these $k_n$s over the ends of the function (being + and - infinity) what we are essentially doing is creating a grid in the $x$ axis with interval $$\Delta k = \frac{\pi}{L}$$ We define function $$F(a) = \int^L_{-L} f(x) e^{ia x}$$ Using this, we can rewrite our complex coefficient $c_n$ to be $$c_n = \frac{1}{2L}F(k_n)$$ Therefore, the Fourier exponential series with $c_n$ becomes: $$f(x) = \sum^{\infty}_{n=-\infty} \frac{F(k_n)}{2L} e^{ik_nx}$$ Notice, however, that $\frac{1}{2L} = \frac{\Delta k}{2\pi}$. Therefore: $$f(x) = \sum^{\infty}_{n=-\infty} \frac{1}{2\pi}F(k_n) e^{ik_nx} \Delta k,$$ which is a much more aesthetic (and for our purposes useful) form. Recall that a Riemann sum approximates the integral of a function $$\int^b_af(x)$$ as the areas of a set of rectangles of equal width under a curve as the number of rectangles tends to infinity (or conversely the width of the rectangles tends to 0). It is written as: $$\sum_{i=1}^{n} f(x_i^*) \Delta x$$

where $\Delta x$ is the width of each subinterval, and $f(x_i^*)$ represents the value of the function at a particular point $x_i^*$ within each subinterval (this can be done in terms of left endpoint, right endpoint and centrepoint sums, being in what way the rectangle is aligned with the curve, however for our purposes we deal only with the left endpoint Riemann sum). The similarity between our summation and the Riemann sum should now become visible - consider what happens as L approaches infinity in our previous summation: $$f(x) = \sum^{\infty}_{n=-\infty} \frac{1}{2\pi}F(k_n) e^{ik_nx} \Delta k.$$ We know that $k_n$ will go to zero. Therefore, $e^{ik_nx}$ will go to 1. Recall the grid image I described earlier - as $\Delta k$ goes to 0, the grid will become infinitely dense. We can take out the $\frac{1}{2\pi}$ to yield $$f(x) = \frac{1}{2\pi}\sum ^{\infty}_{n=-\infty} F(k_n) e^{ik_nx} \Delta k.$$ Now, if $\Delta k$ goes to zero (which it does), the Riemann sum can be made equivalent to the integral by saying that $$f(x) =_{L \to \infty} \frac{1}{2\pi}\int ^{\infty}_{-\infty} F(k)e^{ikx} \, dk$$ for $$F(k) = \int^{\infty}_{-\infty} f(x) e^{ikx} \, dx$$ Which extends the Fourier series across the entire real number line.

Chebyshev Approximation

A method of approximation came to the great Russian mathematician P. L. Chebyshev whilst he was working on the problem of linkages, specifically relating to translating the linear motion of steam engines to the circular motion of wheels. Linkages are devices composed of rigid bodies (called links) connected by joints or pivots. The problem of linkages more generally is that of converting any geometry of motion into any other geometry, which must not necessarily be linear. The problem is still greatly applicable today, examples of its application being in robotics, prosthetics and linkage based engines (which is the problem Chebyshev was working on), and research is still in progress, examples of papers including https://doi.org/10.48550/arXiv.1312.1060 (An algebraic study of linkages with helical joints, Abban et al., 2024) and https://doi.org/10.48550/arXiv.2307.13839 (Geometry of integrable linkages, Perline and Tabachnikov, 2023). Chebyshev is often called the father of Russian mathematics, a claim which is particularly backed up by his founding of the St. Petersburg Mathematical School, kickstarting Russian mathematics.

Pascal’s triangle and Chebyshev polynomials

Pascal’s triangle is formed by beginning with 1 at the top, and summing each number with the number next to it (imagining the triangle is in a space of infinite zeroes) forming the triangular array of the binomial coefficients arising in probability theory, combinatorics, and algebra. The beginning of the triangle is shown above as Figure 1.2. Pascal’s triangle is intrinsically linked to the Fibonacci numbers, a sequence defined inductively as: $$F_0 = 0, \quad F_1 = 1$$ $$F_n = F_{n-1} + F_{n-2}, \quad \text{for } n > 1$$ and with the first few terms being: $$0, 1, 1, 2, 3, 5, \ldots$$

Diagram reconstructed for this web edition; original image not supplied.

Figure 1.2 also shows red lines, travelling in a northeast pattern and each time travelling down one place and to the right two from the previous index of the numbers. The numbers on the red line are then summated, resulting in the Fibonacci sequence. From these lines, we find the coefficients of an interesting set of polynomials: the Fibonacci polynomials. These are polynomials which use the same recursion as the Fibonacci numbers, and add a variable. They are defined as: $$F_n(x) = \begin{cases} 0 & \text{if } n = 0 \\ 1 & \text{if } n = 1 \\ xF_{n-1}(x) + F_{n-2}(x) & \text{if } n \geq 2 \end{cases}$$ We can therefore write out the first few Fibonacci polynomials as: $$\begin{align*} F_0(x) &= 0 \\ F_1(x) &= 1 \\ F_2(x) &= x \\ F_3(x) &= x^2 + 1 \\ F_4(x) &= x^3 + 2x \\ F_5(x) &= x^4 + 3x^2 + 1 \\ F_6(x) &= x^5 + 4x^3 + 3x \\ \end{align*}$$ Notice that the coefficients of the $n^{th}$ Fibonacci polynomial are simply the numbers along the $n^{th}$ red line. We can play around with Pascal’s triangle by changing the signs in every second diagonal. $$\begin{array}{ccccccccccccccc} & & & & & & & 1 & & & & & & & \\ & & & & & & 1 & & -1 & & & & & & \\ & & & & & 1 & & -2 & & 1 & & & & & \\ & & & & 1 & & -3 & & 3 & & -1 & & & & \\ & & & 1 & & -4 & & 6 & & -4 & & 1 & & & \\ & & 1 & & -5 & & 10 & & -10 & & 5 & & -1 & & \\ & 1 & & -6 & & 15 & & -20 & & 15 & & -6 & & 1 & \\ 1 & & -7 & & 21 & & -35 & & 35 & & -21 & & 7 & & -1 \\ \end{array}$$

From this, we can form a new set of polynomials: the Chebyshev polynomials of the second kind. These are formed in a very similar way as we formed the Fibonacci polynomials: $$U_n(x) = \begin{cases} 0 & \text{if } n = 0 \\ 1 & \text{if } n = 1 \\ xU_{n-1}(x) - U_{n-2}(x) & \text{if } n \geq 2 \end{cases}$$ this definition being identical to the one of Fibonacci polynomials, save that we changed the sign from a + to a - in the expression $xU_{n-1}(x) - U_{n-2}(x)$. As well as the Chebyshev polynomials of the second kind, Chebyshev also introduced the Chebyshev polynomials of the first kind, which are the polynomials which uniquely satisfy conditions $${\displaystyle T_{n}(x)={\begin{cases}\cos(n\arccos x)&{\text{ if }}~|x|\leq 1\\\cosh(n\operatorname {arcosh} x)&{\text{ if }}~x\geq 1\\(-1)^{n}\cosh(n\operatorname {arcosh} (-x))&{\text{ if }}~x\leq -1\end{cases}}}$$ The Chebyshev,polynomials of the first kind are also highly useful for approximation, but we will focus here on the second kind. Here the coefficient operate in the same way, being the same diagonal lines on the triangle as formed before (unfortunately LaTeX does not allow me to format the lines on a non-figure). Given this method of computing them, we can consider the first few Chebyshev polynomials and their graph below: $$\begin{align*} U_0(x) &= 1 \\ U_1(x) &= x \\ U_2(x) &= 2x^2 - 1 \\ U_3(x) &= 4x^3 - 3x \\ U_4(x) &= 8x^4 - 8x^2 + 1 \\ \end{align*}$$

This graph is neat (a lot more so than that of the Fibonacci polynomials), and shows one of the fundamental properties of the Chebyshev polynomials, which will be of use to us later: they have very neat and roots.

Using the Chebyshev polynomials for approximation

These Chebyshev polynomials are extremely useful for approximation. This section will concern itself only with the Chebyshev polynomials of the second kind. What is important about them is that they minimize oscillations that can occur at the edges of the interval when using polynomial approximations. This is why, when Chebyshev polynomials are used for approximations, they are commonly only used between the interval [-1, 1]. The reason that the interval [-1, 1] is perfect for the use of Chebyshev polynomials is because Chebyshev polynomials $U_n(x)$ have exactly $n$ real roots evenly (due to the symmmetry of $U_n(x)$ about the origin) spaced out (as is visible in the graph above). Therefore, on the interval [-1, 1], function $f(x)$ can be represented as $$f(x) = \sum^{\infty}_{k=0} a_k U_k(x)$$ Similarly as in the previous section on Fourier analysis, we need to derive $a_k$. We will again, as in the section on Fourier analysis, need to exploit the property of orthogonality. Recall that two functions are orthogonal if $$\int^L_{-L}f(x)g(x)w(x) \, dx= 0$$ and we let $w(x)$ be the weight function (which we have previously just kept as 1). Now consider whether this is true in the case of the Chebyshev polynomials. Consider $$\langle U_m(x)U_n(x)\rangle= \int^1_{-1} U_m(x)U_n(x)$$ The Chebyshev polynomials are constructed in such a way that they are orthogonal with respect to weight function $w(x) = \frac{1}{\sqrt{1-x^2}}$ over the interval $[-1, 1]$. Therefore, we can write that:

$$\int^1_{-1} \frac{U_m(x)U_n(x)}{1 - x^2} = 0, m \neq n$$ And for $m=n$, $$\langle U_m(x)U_n(x)\rangle = \int^1_{-1}U^2_m(x),$$ which yields (with respect to $w(x)$): $$\langle U_m, U_n \rangle = \int^1_{-1} U_m(x)U_n(x) = \begin{cases} 0 & \text{if } m \neq n \\ \pi & \text{if } m = n \neq 0 \\ \pi/2 & \text{if } m = n = 0 \end{cases}$$ Having now demonstrated that the Chebyshev polynomials are orthogonal, we can now levy this property by multiplying our previous expression by another Chebyshev polynomial $U_m$. Therefore $$f(x)U_m(x) = (\sum^{\infty}_{k=0}a_kU_k(x)(x))U_m(x)$$ We can integrate over [-1, 1] so that $$\int^{1}_{-1}f(x)U_m(x)\, dx = \int^1_{-1}(\sum^{\infty}_{k=0}a_kU_k(x)(x))U_m(x) \, dx$$ Since we have discovered the behaviour of the integral of the multiple of two Chebyshev polynomials, we can say that the integral is 0, except when $m = k$, $$\int^1_{-1}f(x)U_k(x)(x)\, dx = a_k \int^1_{-1}U^2_m(x)\,dx$$ Since we know that $$\int^1_{-1} U_k(x)^2(x)= \begin{cases}\pi & \text{if } m = n \neq 0 \\ \pi/2 & \text{if } m = n = 0 \end{cases}$$ we can say that $$\int_{-1}^{1} f(x) U_kx) \, dx = \begin{cases} \pi a_k & \text{if } k \neq 0 \\ \pi/2 a_k & \text{if } k = 0 \end{cases}.$$ Therefore, we divide by the respective constant to yield: $$a_k= \frac{1}{\pi} \int_{-1}^{1} f(x)U_k(x)(x) \, dx, k \neq 0$$ And $$a_0 = \frac{1}{2} \int_{-1}^{1} f(x) \, dx$$ Note that we need to multiply these constants by the weight function $w(x) = \frac{1}{\sqrt{1-x^2}}$ discussed earlier, in order to ensure the integral of the inner product can be correctly calculated with the orthogonality. Therefore: $$a_k = \frac{2}{\pi} \int_{-1}^{1} \frac{f(x) U_k(x)(x)}{\sqrt{1-x^2}} \, dx, \quad \text{for } k \neq 0$$ and $$a_0 = \frac{1}{\pi} \int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} \, dx$$

Approximation of $\pi$

Mathematics is filled with beautiful and important irrational constants. This chapter will concern itself with just $\pi$. $\pi$ is one of the most fundamental constants in all of mathematics (and I think the only mathematical object to which one of the 365 days of the year is dedicated) - it is the ratio formed by the geometric length of a locus of points equidistant with radius $r$ from a centrepoint, and the length $2r$. It is therefore the ratio between the diameter and the circumference of a circle. Hence comes the formula for the circumference $$c = \pi d.$$ This value was first proposed by Archimedes. However, how can we approximate it? How do we know it needs to be approximated and is not just rational?

Hermite’s proof that $\pi$ is irrational

A key question in the consideration of whether or not a constant can be approximated to complete accuracy is the question of whether or not it is irrational. This question was first considered by Johann Heinrich Lambert in the 1760s, who proved that $\pi$ is irrational, that is to say it cannot be described by a proper fraction. Lambert demonstrated this (although I will not go into the particulars here) by showing that this continued fraction expansion holds: $$\tan(x)=\frac{x}{1-\frac{x^2}{3-\frac{x^2}{5-\frac{x^2}{7-\ddots}}}}.$$

Then Lambert proved that if $x$is non-zero and rational, then this expression must be irrational. However, $\tan \frac{\pi}{4} = 1$, therefore $\frac{\pi}{4}$ is irrational, and so, therefore, is $\pi$.

However, in the 19th century Charles Hermite created a much simpler proof, relying only on elementary calculus, such as that taught in the last few years of secondary schools today. Simplicity of ideas used is key in the development of fundamental proofs, and I think that Hermite’s proof demonstrates this beautifully. Hermite’s proof is now the most commonly used and taught proof for $\pi$’s irrationality.

Hermite’s proof operates through use of Reductio ad Absurdem, one of the most powerful techniques of proof by contradiction. It’s power was put perfectly by mathematician GH Hardy in his book A Mathematician’s Apology: "a chess player may offer the sacrifice of a pawn or even a piece, but a mathematician offers the game." (many thanks to NA for gifting me said book.) Let us first assume that $\pi$ can be expressed by a fraction: $$\pi = \frac{u}{v}$$ The next step of the proof is to define a sequence. This sequence is very interesting - its terms are as follows: $$\begin{align*} I_0 &= 2 \\ I_1 &= 4v^2 \\ &\vdots \\ I_n &= \frac{v^{2n}}{n!}\int_0^\pi x^n(\pi - x)^n\sin x \, dx \end{align*}$$

I realise that at first, resulting out of such simple terms $I_0$ and $I_1$ we have such an immensely complicated $n^{th}$ term $I_n$, however the explanation for the formula is remarkably simple, and the resulting consequences are immense. Let us begin by attempting to evaluate the integral $\int^\pi_0 x^n(\pi - x)^n\sin x dx$. Consider the function $f(x, n) = x^n(\pi - x)^n\sin x$ being integrated, defined over the interval [0, $\pi$] with variable positive values for n. The graphs all look relatively similar to this one (where n = 0.5):

Note the fact that the curve begins at 0, comes up to a maximum value and then goes down to 0 again. Let us now consider the maximum possible value for the function, to get a rough idea for the area of the curve. Since $0 \leq x \leq \pi$ we know that the factor $x^n$ in the integral cannot be greater than $\pi^n$, achieved when $x = \pi$. The second term is similarly equal to $\pi^n$, however this is achieved when $x = 0$. Furthermore, $\sin x$ can be at most 1. Therefore, when we take into account the entire function which is being integrated, we can definitively say that: $$f(x, n) \leq \pi^{2n}$$ Therefore, the maximum value for the function must be below the line $y = \pi^{2n}$. From this we can extend that the area $\int^\pi_0 x^n(\pi - x)^n\sin x dx$ must be less than the area formed by the lines $y = \pi^{2n}$ and $x = \pi$, $x = 0$ and $y = 0$. Therefore, $$0 < I_n \leq \frac{(v)^{2n}}{n!}\pi^{2n}\pi$$ Simplification yields that: $$0 < I_n \leq \frac{(v\pi)^{2n}}{n!}\pi$$ We now consider what happens when $n \to \infty$. Let us primarily consider the fraction $\frac{(u\pi)^{2n}}{n!}$, if we let $(u\pi)^2 = x$, we are left with the more general problem - to what value will $\frac{x^n}{n!}$ tend as $n \to \infty$? As we now that factorial growth is greater than exponential growth, it is demonstrated that it will, in fact converge to 0. Therefore, it follows that sequence $$I_0, I_1, I_2, ... I_N \to 0.$$ This is a fundamental step in the proof that $\pi$ is irrational. Recall that we begun with the assumption that $\pi$ can be described as a fraction $\frac{u}{v}$, where $u$ and $v$ are positive integers. We have now demonstrated that the sequence $I_0, I_1, I_2, ... I_N$ converges to 0, if we can prove that this sequence is also made up of only positive integers based on the assumption that $u$ and $v$ are positive integers, we will have an absurdity - positive integers cannot converge to 0! Therefore we will have proved that $u$ and $v$ cannot be positive integers, proving that $\pi$ is irrational. In order to establish the fact that all integers in $I$ are positive integers, we need to create a recurrence relation for $I_n$, based on $I_{n-1}$ and $I_{n-2}$.

The proof that all values of the sequence lies in a way of rewriting this sequence as an inductive sequence. Earlier, I stated the first to values of the sequence I as $$I_0 = 2$$ and $$I_1 = 4v^2$$

Let me quickly demonstrate that these are true. The first is just the result of the simple and easy to evaluate integral $$I_0 = \int^\pi_0 \sin x dx = 2,$$ The second term is slightly more difficult to derive, it comes from the integral $$I_1 = v^2 \int^\pi_0 x(\pi - x) \sin x dx.$$ This simplifies to $$I_1 = v^2 \pi \int^\pi_0 x \sin x dx - v^2\int^\pi_0 x^2 \sin x dx$$ This simplified integral simplifies to two much simpler, standard integrals for integration by parts, which simplifies the integral further to $$I_1 = v^2\pi^2 - v^2\pi^2 + 4v^2,$$ which simplifies to $$I_1 = 4v^2.$$

Now we have fully established the starting values of the sequence, let us again use integration by parts to establish a recurrence relation. Let us evaluate the integral on a specific case in the sequence, namely that for $I_5$. Consider the integral $$\int_0^\pi x^5(\pi - x)^5\sin xdx$$ We can simplify this to: $$\int_0^\pi 5(\pi - 2x)x^4(\pi - x)^4\cos x \, dx$$ Note the interesting property that the formula looks very similar to what we established as the nth term, only the exponents are down by one, sine is replaced with cosine and we have some extra terms at the start. Let us now integrate by parts again, resulting in $$\int_0^\pi (10x^4(\pi - x)^4 - 20(\pi^2 - 4\pi x + 4x^2)x^3(\pi - x)^3)\sin x\,dx$$ Some basic algebra simplifies this to $$\int_0^\pi -20 \pi^2x^3(\pi - x)^3 + 90x^4(\pi - x)^4$$ This presents important progress, because the above expression can be split into two separate integrals, resulting in $$-20 \pi^2\int^\pi_0 x^3 (\pi - x)^3 \sin x \, dx + 90\int^\pi_0 x^4(\pi - x)^4 \sin x \, dx$$

This looks quite neatly simplified so far, so let’s now input it into our general formula for $I_5$: $$I_5 = -\frac{v^{10}}{5!}20 \pi^2\int^\pi_0 x^3 (\pi - x)^3 \sin x \, dx + \frac{v^{10}}{5!}90\int^\pi_0 x^4(\pi - x)^4 \sin x \, dx$$ This gives us one final step of simplification, resulting in $$I_5 = -u^4 \pi^2 \frac{v^6}{3!}\int^\pi_0 x^3(\pi - x)^3\sin x \, dx + 2v^2 9 \frac{v^8}{4!}\int^\pi_0 x^4 (\pi - x)^4 \sin x \, dx.$$ This results in our proof. We can now simply replace $\frac{v^6}{3!}\int^\pi_0 x^3(\pi - x)^3\sin x$ with $I_3$, and $\frac{v^8}{4!}\int^\pi_0 x^4 (\pi - x)^4 \sin x \, dx$ with $I_4$. Off the back of this, we can build an inductive sequence, $$I_5 = 18v^2 I_4 - v^4\pi^2I_3.$$ We can now make use of our assumption that pi is rational: $$I_5 = 18v^2 I_4 - v^2u^2I_3$$ Also note that $18 = 2(2 \cdot 5 \cdot - 1)$. If we wanted to generalise this, we could write that $$I_n = 2v^2(2n - 1)I_{n-1} - v^2u^2I_{n-2}$$ Note that this does not prove the relation for all $n$. In order to prove this for all $n$, we can begin by considering the sequences of real functions: $$\begin{aligned}A_{0}(x)&=\sin(x),&&A_{n+1}(x)=\int _{0}^{x}yA_{n}(y)\,dy\\[4pt]U_{0}(x)&={\frac {\sin(x)}{x}},&&U_{n+1}(x)=-{\frac {U_{n}'(x)}{x}}\end{aligned}$$

We can inductively prove that these functions can be represented as infinite series, such that $${\displaystyle {\begin{aligned}A_{n}(x)&={\frac {x^{2n+1}}{(2n+1)!!}}-{\frac {x^{2n+3}}{2\times (2n+3)!!}}+{\frac {x^{2n+5}}{2\times 4\times (2n+5)!!}}\mp \cdots \\[4pt]U_{n}(x)&={\frac {1}{(2n+1)!!}}-{\frac {x^{2}}{2\times (2n+3)!!}}+{\frac {x^{4}}{2\times 4\times (2n+5)!!}}\mp \cdots \end{aligned}}}$$ and therefore that $$U_{n}(x)={\frac {A_{n}(x)}{x^{2n+1}}}.\,$$ From this, we can say that (using the quotient rule for differentiation) $${\displaystyle {\begin{aligned}{\frac {A_{n+1}(x)}{x^{2n+3}}}&=U_{n+1}(x)=-{\frac {U_{n}'(x)}{x}}=-{\frac {1}{x}}{\frac {\mathrm {d} }{\mathrm {d} x}}\left({\frac {A_{n}(x)}{x^{2n+1}}}\right)\\[6pt]&=-{\frac {1}{x}}\left({\frac {A_{n}'(x)\cdot x^{2n+1}-(2n+1)x^{2n}A_{n}(x)}{x^{2(2n+1)}}}\right)\\[6pt]&={\frac {(2n+1)A_{n}(x)-xA_{n}'(x)}{x^{2n+3}}}\end{aligned}}}$$ This is equivalent to saying that $${\displaystyle A_{n+1}(x)=(2n+1)A_{n}(x)-x^{2}A_{n-1}(x).\,}$$ Notice an inherent similarity between this recurrence relation and the one which was demonstrated for $I_5$ - the final step in the proof is the idea that we can demonstrate that $A_n$ and $I_n$ are essentially the same function.