Skip to main content

Section 26.3 Cosine interpolation of non-periodic functions

This is an optional section, feel free to skip it.

In Chapter 25 we worked with periodic functions, approximating them by trigonometric polynomials. Can we use trigonometric polynomials to approximate non-periodic functions? For example, suppose we solved the initial-value problem

\begin{equation} y'(t) =\sin(y^2)+t,\quad y(0)=1\label{eq-diffeq-for-intepolation}\tag{26.3.1} \end{equation}

with

[t, y] = ode45(@(t, y) sin(y^2) + t, 0:0.1:3, 1);

Note that in contrast to Chapter 20 the second argument is given not as an interval (two-element vector) [0 3] but as a set of t-values at which the solution needs to be found. This ensures that we get y-values for these specific t-values. But now that we have them, how can we get a simple approximate solution?

  • Polynomial interpolation? Hopeless, a degree 30 polynomial through 31 equally spaced points is not even worth trying.
  • Cubic spline interpolation? It would make a nice plot but in terms of a formula, the spline would have 120 coefficients (4 for each of 30 subintervals). This is hardly a simple solution.
  • Trigonometric polynomials? But they are periodic, and our function has different values at 0 and 3.

Compare the behavior of cosines \(\cos kx\) and sines \(\sin kx\) on the interval \([0, \pi]\text{.}\) The cosines have different values at \(0, \pi\) while the sines are always 0 at both ends. This makes cosines, specifically in the form \(\cos (\pi k x/L)\text{,}\) preferable for approximating a non-periodic function on \([0, L]\text{.}\) The Discrete Cosine Transform is a tool used to find the coefficients of a given function in the cosine basis. It is somewhat technical (there are 4 different types of the transform), but the essence of the method can be stated simply: apply Fourier transform to the even periodic extension of the given function.

The discrete version of even periodic extension is as follows: given function values \(y_1, \dots, y_n\text{,}\) we extend them to \(y_1, \dots, y_n, y_{n-1}, \dots, y_2\) (for the total of \(2n-2\) values). In Matlab terms, the extended vector is

ye = [y y(end-1:-1:2)];

Why not simply [y y(end:-1:1)]? Repeating the endpoint values \(y_1, y_n\) would distort the picture of this function. They should appear once while the rest appear twice; this can be viewed as yet another form of “trapezoidal adjustment”.

Having found the Fourier coefficients c = fft(ye)/numel(ye) we observe they are all real (any imaginary component is due to round-off errors). This makes it possible to “fold” the sum of complex exponentials into a cosine sum: for \(k=1, \dots, n-2\) the pair of terms

\begin{equation*} c_k \exp(i k x ) + c_{2n-k} \exp(i (2n-k) x ) \end{equation*}

is aliased to

\begin{equation*} c_k \exp(i k x ) + c_{2n-k} \exp(i (-k) x ) = 2c_k \cos (kx) \end{equation*}

(here \(x\) is a multiple of \(\pi/n\text{.}\)) The terms \(c_0\) and \(c_{n-1}\) do not have a counterpart among the coefficients \(c_0, \dots, c_{2n-2}\text{,}\) so the cosine sum ends up being

\begin{equation} c_0 + 2\sum_{k=1}^{n-2} c_k \cos (\pi kx/L) + c_{n-1} \cos (\pi (n-1) x/L)\label{eq-cosine-sum}\tag{26.3.2} \end{equation}

And of course, Matlab indexing of the coefficients begins with 1, forcing another adjustment: see the following example.

Apply the cosine interpolation to the numerical solution of ODE (26.3.1).

Solution

Since ode45 returns y as a column, it is transposed in the code below. Then we extend this vector to ye, take its Fourier Transform c, and fold it into cosine coefficients cf. Then plot the resulting sum of cosines.

[t, y] = ode45(@(t, y) sin(y^2) + t, 0:0.1:3, 1);
y = y';

n = numel(y);
ye = [y y(end-1:-1:2)];
c = real(fft(ye))/numel(ye);
cf = [c(1) 2*c(2:n-1) c(n)];

L = 3;
x = linspace(0, L, 1000);
k = 0:n-1;
g = cf*cos(k'*x*pi/L);
plot(x, g, t, y, 'r*')

The cosine sum in Example 26.3.1 has 31 terms, as many as there are data points. But many of them are very close to zero, as plot(cf) shows. For example, we could keep the first m=8 of them and still have a good approximation for this function.

m = 8;
k = 0:m-1;
g = cf(1:m)*cos(k'*x*pi/L);

This is essentially a one-dimensional version of image compression discussed in Section 26.2.