Trigonometric interpolation - Fundamentals of Numerical Computation

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)

Julia

MATLAB

Python

Trigonometric interpolation

triginterp.jl

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

Construct the trigonometric interpolant for the points defined by
vectors `t` and `y`.
"""
function triginterp(t, y)
    N = length(t)

    τ(x) =
        if x == 0
            return 1.0
        else
            denom = isodd(N) ? N * sin(π * x / 2) : N * tan(π * x / 2)
            return sin(N * π * x / 2) / denom
        end

    return function (x)
        return sum(y[k] * τ(x - t[k]) for k in eachindex(y))
    end
end

Trigonometric interpolation

triginterp.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function p = triginterp(t, y)
% TRIGINTERP Trigonometric interpolation.
% Input:
%   t   equispaced interpolation nodes (vector, length N)
%   y   interpolation values (vector, length N)
% Output:
%   p   trigonometric interpolant (function)

N = length(t);
p = @value;

    function f = value(x)
        f = zeros(size(x));
        for k = 1:N
            f = f + y(k) * trigcardinal(x - t(k));
        end
    end    % value function
        
    function tau = trigcardinal(x)
        if rem(N,2)==1   % odd
            tau = sin(N * pi * x / 2) ./ (N * sin(pi * x / 2));
        else             % even
            tau = sin(N * pi * x / 2) ./ (N * tan(pi * x / 2));
        end
        tau(isnan(tau)) = 1;    % fix any divisions by zero (L'Hopital's Rule)
    end    % trigcardinal function
end

Trigonometric interpolation

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)

Example 9.5.1 (Trigonometric interpolation)

Julia

MATLAB

Python

Example 9.5.1

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.

using Plots
N = 7
n = (N - 1) ÷ 2
t = 2 * (-n:n) / N
y = zeros(N)
y[n+1] = 1

p = FNC.triginterp(t, y);
plot(p, -1, 1)

scatter!(t, y, color=:black, 
    xaxis=(L"x"),  yaxis=(L"\tau(x)"),
    title="Trig cardinal function, N=$N")

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

f(x) = exp( sinpi(x) - 2*cospi(x) )
y = f.(t)
p = FNC.triginterp(t, y)

plot(f, -1, 1, label="function",
    xaxis=(L"x"),  yaxis=(L"p(x)"),
    title="Trig interpolation, N=$N", legend=:top)
scatter!(t, y, m=:o, color=:black, label="nodes")
plot!(p, -1, 1, label="interpolant")

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.

Source

N = 2:2:60
err = zeros(size(N))
x = range(-1, 1, 2501)  # for measuring error
for (k,N) in enumerate(N)
    n = (N-1) / 2;   t = 2*(-n:n) / N;
    p = FNC.triginterp(t, f.(t))
    err[k] = norm(f.(x) - p.(x), Inf)
end

plot(N, err, m=:o,
    xaxis=(L"N"),  yaxis=(:log10, "max error"),
    title="Convergence of trig interpolation")

Example 9.5.1

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.

N = 7;  n = (N-1) / 2;
t = 2 * (-n:n)' / N;
y = zeros(N, 1);  y(n+1) = 1;
clf,  scatter(t, y, 'k'),  hold on

p = triginterp(t, y);
fplot(p, [-1, 1])
xlabel('x'),  ylabel('p(x)')   
title('Trig cardinal function')

Unrecognized function or variable 'triginterp'.

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

clf
f = @(x) exp( sin(pi*x) - 2 * cos(pi*x) );
fplot(f, [-1, 1], displayname="periodic function"),  hold on
fplot(triginterp(t, f(t)), [-1, 1], displayname="trig interpolant")
y = f(t);  scatter(t, f(t), 'k')
xlabel('x'),  ylabel('f(x)')   
title('Trig interpolation');  legend()

Unrecognized function or variable 'triginterp'.

Source

N = 2:2:60;
err = zeros(size(N));
x = linspace(-1, 1, 1601)';  % for measuring error
for k = 1:length(N)
    n = (N(k) - 1) / 2;
    t = 2 * (-n:n)' / N(k);
    p = triginterp(t, f(t));
    err(k) = norm(f(x) - p(x), Inf);
end
clf,  semilogy(N, err, 'o-')
axis tight, title('Convergence of trig interpolation')   
xlabel('N'),  ylabel('max error')

Unrecognized function or variable 'triginterp'.

Example 9.5.1

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.

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");

Source

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{a_0}{2}, & 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)

Julia

MATLAB

Python

Example 9.5.2

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

f(x) = 3 * cospi(2x) - cispi(x)    # cispi(x) := exp(1im * π * x)

f (generic function with 1 method)

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

n = 4;  N = 2n+1;
t = [ 2j / N for j in 0:N-1 ]      # nodes in [0,2)
y = f.(t);

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

using FFTW
c = fft(y) / N
freq = [0:n; -n:-1]
@pt :header=["k", "coefficient"] [freq round.(c, sigdigits=5)]

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) )$ .

Source

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

freq = [0:n; -n:-1]
scatter(freq, abs.(c);
    xaxis=(L"k", [-n, n]),  yaxis=(L"|c_k|", :log10), 
    title="Fourier coefficients",  legend=:none)

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.

Example 9.5.2

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

f = @(x) 3 * cos(2*pi * x) - exp(1i*pi * x);

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

n = 4;
N = 2*n + 1;
t = 2 * (0:N-1)' / N;      % nodes in $[0,2)$
y = f(t);

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

c = fft(y) / N;
freq = [0:n, -n:-1]';
format short
table(freq, c, variableNames=["k", "coefficient"])

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

Source

f = @(x) exp( sin(pi*x) );    % content at all frequencies
n = 9;  N = 2*n + 1;
t = 2 * (0:N-1)' / N;         % nodes in $[0,2)$
y = f(t);
c = fft(y) / N;
freq = [0:n, -n:-1]';

clf
semilogy(freq, abs(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.

Example 9.5.2

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

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.

Source

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=-\epsilon_\text{mach}$ . 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.

Preface

Orthogonal polynomials

Preface

Spectrally accurate integration