from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

9.5.Trigonometric interpolation¶

Up to this point, all of our global approximating functions have been polynomials. While they are versatile and easy to work with, they are not always the best choice.

Suppose we want to approximate a function $f$ that is periodic, with one period represented by the standard interval $[-1,1]$ . Mathematically, periodicity means that $f(x+2)=f(x)$ for all real $x$ . We could use polynomials to interpolate or project $f$ . However, it seems more reasonable to replace polynomials by functions that are also periodic.

It turns out that trigonometric interpolation allows us to return to equally spaced nodes without any problems. We therefore define $N=2n+1$ equally spaced nodes inside the interval $[-1,1]$ by

t_k = \frac{2k}{N}, \quad k=-n,\ldots,n.

(9.5.2)

The formulas in this section require some minor but important adjustments if $N$ is even instead. We have modified our standard indexing scheme here to make the symmetry within $[-1,1]$ about $x=0$ more transparent. Note that the endpoints $\pm 1$ are not among the nodes.

As usual, we have sample values $y_{-n},\ldots,y_n$ , perhaps representing values of a function $f(x)$ at the nodes. We also now assume that the sample values can be extended periodically forever in both directions, so that $y_{k+mN}=y_k$ for any integer $m$ .

9.5.1Cardinal functions¶

We can explicitly state the cardinal function basis for equispaced trigonometric interpolation. It starts with

\tau(x) = \frac{2}{N} \left( \frac{1}{2} + \cos \pi x + \cos 2\pi x + \cdots + \cos n\pi x\right) = \frac{\sin(N\pi x/2)}{N\sin(\pi x/2)}.

(9.5.3)

You can directly check the following facts. (See Exercise 9.5.3.)

Because the functions $\tau_{-n},\ldots,\tau_n$ form a cardinal basis, the coefficients of the interpolant are just the sampled function values, i.e., the interpolant of points $(t_k,y_k)$ is

p(x) = \sum_{k=-n}^n y_k \tau_k(x).

(9.5.4)

The convergence of a trigonometric interpolant is spectral, i.e., exponential as a function of $N$ in the max-norm.

9.5.2Implementation¶

Function 9.5.1 is an implementation of trigonometric interpolation based on (9.5.4). The function accepts an $N$ -vector of equally spaced nodes. Although we have not given the formulas above, the case of even $N$ is included in the code.

Algorithm 9.5.1 (triginterp)

triginterp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def triginterp(t, y):
    """
        triginterp(t, y)

    Return trigonometric interpolant for points defined by vectors t and y.
    """
    N = len(t)

    def trigcardinal(x):
        if x == 0:
            tau = 1.0
        elif np.mod(N, 2) == 1:  # odd
            tau = np.sin(N * np.pi * x / 2) / (N * np.sin(np.pi * x / 2))
        else:  # even
            tau = np.sin(N * np.pi * x / 2) / (N * np.tan(np.pi * x / 2))
        return tau

    def p(x):
        return np.sum([y[k] * trigcardinal(x - t[k]) for k in range(N)])

    return np.vectorize(p)

About the code

The construct on line 13 is known as a ternary operator. It is a shorthand for an if–else statement, giving two alternative results for the true/false cases. Line 19 uses eachindex(y), which generalizes 1:length(y) to cases where a vector might have a more exotic form of indexing.

Example 9.5.1 (Trigonometric interpolation)

We will get a cardinal function without using an explicit formula, just by passing data that is 1 at one node and 0 at the others.

Tip

The operator ÷, typed as \div then Tab, returns the quotient without remainder of two integers.

N = 7
n = int((N - 1) / 2)
t = 2 * arange(-n, n + 1) / N
y = zeros(N)
y[n] = 1

p = FNC.triginterp(t, y)
x = linspace(-1, 1, 600)
plot(x, p(x))
plot(t, y, "ko")

xlabel("x"),  ylabel("tau(x)")
title("Trig cardinal function");

Here is a 2-periodic function and one of its interpolants.

f = lambda x: exp(sin(pi * x) - 2 * cos(pi * x))

plot(x, f(x), label="periodic function")
y = f(t)

p = FNC.triginterp(t, y)
plot(x, p(x), label="trig interpolant")
plot(t, y, "ko", label="nodes")

xlabel("$x$"),  ylabel("$p(x)$")
legend(),  title("Trig interpolation");

The convergence of the interpolant is spectral. We let $N$ go needlessly large here in order to demonstrate that unlike polynomials, trigonometric interpolation is stable on equally spaced nodes. Note that when $N$ is even, the value of $n$ is not an integer but works fine for defining the nodes.

N = arange(2, 62, 2)
err = zeros(N.size)

x = linspace(-1, 1, 1601)    # for measuring error
for k in range(N.size):
    n = (N[k] - 1) / 2
    t = 2 * arange(-n, n + 1) / N[k]
    p = FNC.triginterp(t, f(t))
    err[k] = max(abs(f(x) - p(x)))

semilogy(N, err, "-o")
xlabel("N"),  ylabel("max error")
title("Convergence of trig interpolation");

9.5.3Fast Fourier transform¶

Although the cardinal form of the interpolant is useful and stable, there is a fundamental alternative. It begins with an equivalent complex form of the trigonometric interpolant (9.5.1),

p(x) = \sum_{k=-n}^n c_k e^{ik\pi x}.

(9.5.5)

The connection is made through Euler’s formula,

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

(9.5.6)

and the resultant identities

\cos \theta = \frac{e^{i \theta}+e^{-i\theta}}{2}, \qquad \sin \theta = \frac{e^{i \theta}-e^{-i\theta}}{2i}.

(9.5.7)

Specifically, we have

c_k = \begin{cases} \frac{1}{2} a_0, & k=0, \\[1mm] \frac{1}{2}(a_k - i b_k), & k> 0, \\[1mm] \overline{c_{-k}}, & k < 0. \end{cases}

(9.5.8)

While working with an all-real formulation seems natural when the data are real, the complex-valued version leads to more elegant formulas and is standard.

The $N=2n+1$ coefficients $c_k$ are determined by interpolation nodes at the $N$ nodes within $[-1,1]$ . By evaluating the complex exponential functions at these nodes, we get the $N\times N$ linear system

\mathbf{F}\mathbf{c} = \mathbf{y}, \qquad \mathbf{F} = \bigl[ e^{\,is\pi t_r} \bigr]_{\, r=-n,\ldots,n,\, s=-n,\ldots,n,}

(9.5.9)

to be solved for the coefficients. Up to a scalar factor, the matrix $\mathbf{F}$ is unitary, which implies that the system can be solved in $O(N^2)$ operations simply by a matrix-vector multiplication.

However, one of the most important (though not entirely original) algorithmic observations of the 20th century was that the linear system can be solved in just $O(N\log N)$ operations by an algorithm now known as the fast Fourier transform, or FFT.

The FFTW package provides a function fft to perform this transform, but its conventions are a little different from ours. Instead of nodes in $(-1,1)$ , it expects the nodes to be defined in $[0,2)$ , and it returns the trig polynomial coefficients in the order

\begin{bmatrix} c_0, & c_1, & \cdots & c_n, & c_{-n}, & \cdots & c_{-1} \end{bmatrix}.

(9.5.10)

Example 9.5.2 (FFT)

This function has frequency content at $2\pi$ , $-2\pi$ , and $\pi$ .

f = lambda x: 3 * cos(2 * pi * x) - exp(1j * pi * x)

To use fft, we set up nodes in the interval $[0,2)$ .

n = 4
N = 2 * n + 1
t = 2 * arange(0, N) / N    # nodes in [0,2)
y = f(t)

We perform Fourier analysis using fft and then examine the resulting coefficients.

from scipy.fftpack import fft, ifft, fftshift
c = fft(y) / N
freq = hstack([arange(n+1), arange(-n,0)])
results = PrettyTable()
results.add_column("freq", freq) 
results.add_column("coefficient", c)
results

Note that $1.5 e^{2i\pi x}+1.5 e^{-2i\pi x} = 3 \cos(2\pi x)$ , so this result is sensible.

Fourier’s greatest contribution to mathematics was to point out that every periodic function is just a combination of frequencies—infinitely many of them in general, but truncated for computational use. Here we look at the magnitudes of the coefficients for $f(x) = \exp( \sin(\pi x) )$ .

f = lambda x: exp(sin(pi * x))    # content at all frequencies
n = 9;  N = 2*n + 1;
t = 2 * arange(0, N) / N    # nodes in [0,2)
c = fft(f(t)) / N

semilogy(range(-n, n+1), abs(fftshift(c)), "o")
xlabel("$k$"),  ylabel("$|c_k|$")
title("Fourier coefficients");

The Fourier coefficients of smooth functions decay exponentially in magnitude as a function of the frequency. This decay rate is determines the convergence of the interpolation error.

The theoretical and computational aspects of Fourier analysis are vast and far-reaching. We have given only the briefest of introductions.

9.5.4Exercises¶

Exercise 9.5.5

⌨ Spectral convergence is predicated on having infinitely many continuous derivatives. At the other extreme is a function with a jump discontinuity. Trigonometric interpolation across a jump leads to a lack of convergence altogether, a fact famously known as the Gibbs phenomenon.

(a) Define f(x) = sign(x+eps()). This function jumps from -1 to 1 at $x=-\macheps$ . Plot the function over $-0.05\le x \le 0.15$ .

(b) Let $n=30$ and $N=2n+1$ . Using Function 9.5.1, add a plot of the trigonometric interpolant to $f$ to the graph from part (a).

(c) Repeat part (b) for $n=80$ and $n=180$ .

(d) You should see that the interpolants overshoot and oscillate near the step. The widths of the overshoots decrease with $n$ but the heights approach a limiting value. By zooming in to the graph, find the height of the overshoot to two decimal places.