Chapter 7
Approximation and interpolation
Context: We now apply the idea of basis functions to approximate functions. To do this, we return to the concept of the residual. The goal of function approximation is that the approximated function minimizes the residual. Building on these ideas, we will then discuss the approximation of differential equations in the next chapter.
7.1 Residual
In the previous section, we described how a series expansion can be constructed using basis functions. A typical series expansion contains a finite number of elements \(N\) and has the form \begin {equation} f_N(x) = \sum _{n=1}^N a_n \varphi _n(x), \end {equation} where the \(\varphi _n(x)\) are the basis functions introduced in the previous chapter.
We now want to approach the question of how we can approximate an arbitrary function \(f(x)\) via such a basis function expansion. We define the residual \begin {equation} R(x) = f_N(x) - f(x), \end {equation} which vanishes at every point \(x\) if \(f_N(x)\equiv f(x)\). For an approximation we want to “minimize” this residual. (Minimizing in this context means to bring it as close to zero as possible.) We are looking for the coefficients \(a_n\) of the series, which approximate the function \(f(x)\) in the sense of minimizing the residual.
At this point, it should be noted that the basis functions must be defined on the same support as the target function \(f(x)\). For the approximation of a periodic function \(f(x)\) we need a periodic basis.
7.2 Collocation
The first minimization strategy introduced here is collocation. This method requires that the residual disappears at selected collocation points \(y_n\), \begin {equation} R(y_n) = 0 \quad \text {or}\quad f_N(y_n) = f(y_n). \end {equation} The number of collocation points must correspond to the number of coefficients in the series expansion. The choice of ideal collocation points \(y_n\) itself is non-trivial, and we will only discuss specific cases here.
As a first example, we discuss an expansion into \(N\) finite elements. As collocation points we choose the nodal points of the basis, \(y_n=x_n\). At these sampling points, only one of the basis functions is non-zero, \(\varphi _n(y_n)=1\) and \(\varphi _n(y_k)=0\) if \(n\not =k\). This means that the condition \begin {equation} R(y_n) = 0 \end {equation} trivially leads to \begin {equation} a_n = f(y_n). \end {equation} The coefficients \(a_n\) are therefore the function values at the collocation points. The approximation is a piece-wise linear function between the function values of \(f(x)\).
As a second example, we discuss a Fourier series with corresponding \(N\) Fourier basis functions, \begin {equation} \varphi _n(x) = \exp \left ( i q_n x \right ). \end {equation} In the context of a collocation method, we require that the residual vanishes on \(N\) equidistant points, \(R(y_n)=0\) with \begin {equation} y_n = n L / N, \end {equation} where \(L/N\) is the grid spacing. The collocation condition is \begin {equation} \sum _{k=-(N-1)/2}^{(N-1)/2} a_k \exp \left (i q_k y_n\right ) = \sum _{k=-(N-1)/2}^{(N-1)/2} a_k \exp \left (i 2\pi \frac {k n}{N}\right ) = f(y_n). \label {eq:fourier-collocation} \end {equation} Equations \eqref{eq:fourier-collocation} can now be solved for \(a_k\). We use the fact that for equidistant collocation points the Fourier matrix \begin {equation} W_{kn}=\exp (i2\pi kn/N)=\left [\exp (i2\pi /N)\right ]^{kn} \label {eq:dft-matrix} \end {equation} is unitary (except for a constant factor), i.e. its inverse is given by the adjoint: \begin {equation} \sum _{n=0}^{N-1} W_{kn} W_{nl}^* = \sum _{n=0}^{N-1} \left [\exp (i2\pi /N)\right ]^{n(k-l)} = N\delta _{kl} \end {equation} We can therefore multiply Eq. \eqref{eq:fourier-collocation} by \(W_{nl}^*\) and sum over \(n\). This results in \begin {equation} \sum _{n} \sum _{k} W_{kn} W_{nl}^* a_k = \sum _{k} N a_k \delta _{kl} = N a_l. \end {equation} This means that the coefficients can be expressed as \begin {equation} a_l = \frac {1}{N} \sum _{n=0}^{N} f\left (\frac {nL}{N}\right ) \exp \left (-i 2\pi \frac {l n}{N}\right ) = \frac {1}{N} \sum _{n=0}^{N} f(y_n) \exp \left (-i q_l y_n\right ) \label {eq:inverse-discrete-fourier-transform} \end {equation} for \(-(N-1)\leq l \leq N-1\). This is the discrete Fourier transform (DFT) of the function \(f(y_n)\) discretized on the collocation points.
As a simple example, we show the approximation of the example function \(f(x)=\sin (2\pi x)^3 + \cos (6\pi (x^2-1/2))\) using the Fourier basis and finite elements. Figure 7.2 shows this approximation for \(2N+1=5\) and \(2N+1=11\) basis functions with equidistant collocation points.
The figure shows that all approximations run exactly through the collocation points, as required by the collocation condition. The two approaches interpolate differently between the collocation points. Finite elements lead to a linear interpolation between the points. The Fourier basis is more complicated. The curve between the collocation points is called Fourier interpolation.
7.3 Weighted residuals
We would now like to generalize the collocation method. To do this, we introduce the concept of a test function. Instead of requiring that the residual vanishes at individual points, we require that the scalar product \begin {equation} (v, R) = 0 \label {eq:test-function} \end {equation} with some function \(v(x)\) disappears. If Eq. \eqref{eq:test-function} vanishes for any test function \(v(x)\), then the “weak” formulation Eq. \eqref{eq:test-function} is identical to the strong formulation \(R(x)=0\). Equation \eqref{eq:test-function} is called a “weak” formulation because the condition is only fulfilled in the integral sense. In particular, it is shown later that this weak formulation leads to a weak solution, which cannot satisfy the original (strong) PDGL at every point. The condition \eqref{eq:test-function} is often called a weighted residual, because the residual is weighted by the test function.
A special set of test functions leads directly to the collocation method. We choose the set of \(N\) test functions \begin {equation} v_n(x) = \delta (x-y_n) \label {eq:colloctest} \end {equation} where \(\delta (x)\) is the Dirac \(\delta \) function and \(y_n\) the collocation points. The condition \((v_n,R)=0\) for all \(n\in [0,N-1]\) leads directly to the collocation condition \(R(y_x)=0\).
Note: The Dirac \(\delta \) function should be familiar from lectures on signal processing. The most important property of this function is the filter property, \begin {equation} \int _{-\infty }^{\infty } \dif x\, f(x) \delta (x-x_0) = f(x_0), \end {equation} i.e. the integral over the product of the \(\delta \) function gives the function value at which the argument of the \(\delta \) function disappears. All other properties follow from this, e.g. \begin {equation} \int \dif x\, \delta (x) = \Theta (x), \end {equation} where \(\theta (x)\) is the (Heaviside) step function.
7.4 Galerkin method
The Galerkin method is based on the idea of using the basis functions \(\varphi _n\) of the series expansion as test functions. This leads to the \(N\) conditions \begin {equation} (\varphi _n, R) = 0, \label {eq:galerkinortho} \end {equation} which can be written as \begin {equation} (\varphi _n, f_N) = (\varphi _n, f). \end {equation} For an orthogonal set of basis functions, this yields \begin {equation} a_n = \frac {(\varphi _n, f)}{(\varphi _n, \varphi _n)}. \end {equation} This equation has already been discussed in section ??. For a non-orthogonal basis set, e.g. the basis of the finite elements, the Galerkin condition yields a system of linear equations, \begin {equation} \sum _{m=1}^N (\varphi _n,\varphi _m) a_m = (\varphi _n, f), \label {eq:galerkin-coefficients} \end {equation} where the matrix \(M_{nm}=(\varphi _n,\varphi _m)\) is sparse for the finite elements.
Let us now return to our example function \(f(x)=\sin (2\pi x)^3 + \cos (6\pi (x^2-1/2))\). Figure 7.2 shows the approximation of this function with Fourier and finite element basis sets and the Galerkin method. There are no collocation points and the approximation using finite elements does not exactly match the function to be approximated at the interpolation points. The function is only approximated in the integral sense.
Note: The Galerkin condition (see also Eq. \eqref{eq:galerkinortho}) \begin {equation} (\varphi _n, R) = 0, \end {equation} means that the residual is orthogonal to all basis functions. In other words, the residual can only contain contributions to the function that cannot be described with the given basis set. This implies that we can systematically improve our solution by extending the basis set.
7.5 Least squares
An alternative approach to approximation is to minimize the square of the residual, \((R, R)\), also knows as a least squares approach. For a general series expansion with \(N\) basis functions, we obtain \begin {equation} \begin {split} (R, R) &= (f, f) + (f_N, f_N) - (f_N, f) - (f, f_N) \\ &= (f, f) + \sum _{n=1}^N \sum _{m=1}^N a_n^* a_m (\varphi _n, \varphi _m) - \sum _{n=1}^N a_n^* (\varphi _n, f) - \sum _{n=1}^N a_n (f, \varphi _n). \end {split} \end {equation} This error square is minimized if \begin {equation} \frac {\partial (R,R)}{\partial a_k} = \sum _{n=1}^N a_n^* (\varphi _n, \varphi _k) - (f, \varphi _k) = 0 \end {equation} and \begin {equation} \frac {\partial (R,R)}{\partial a^*_k} = \sum _{n=1}^N a_n (\varphi _k, \varphi _n) - (\varphi _k, f) = 0. \end {equation} This expression is identical to Eq. \eqref{eq:galerkin-coefficients} of the Galerkin method.