Trigonometric interpolation - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *

A smooth function $u : [0,2\pi] \to \re$ can be written in terms of its Fourier series

u(x) = Su(x) = \sum_{k=-\infty}^\infty \hat u_k \phi_k(x), \qquad \phi_k(x) = \ee^{\ii k x}

(1)

where

\hat u_k = \frac{1}{2\pi} \int_0^{2\pi} u(x) \ee^{-\ii k x} \ud x

(2)

is the continuous Fourier transform. Recall that $\{ \phi_k \}$ is an orthogonal family of functions.

We can truncate the series $Su$ to obtain approximations. But this is not practically useful since we cannot compute the coefficients $\hat u_k$ ; either because we may not know $u(x)$ for all $x$ and also because the integrals cannot be computed exactly.

An approximation method which makes use of the function values at a discrete set of points only, will be more useful in computations.

6.12.1Interpolation problem¶

For any even integer $N>0$ , consider the set of points in $[0,2\pi]$

x_j = \frac{2\pi j}{N}, \qquad j=0,1,\ldots,N-1

(3)

which are uniformly spaced with spacing

h = \frac{2\pi}{N}

(4)

N = 9
h = 2*pi/N
x = h * arange(0,N)
figure(figsize=(8,2))
plot(x,0*x,'-o')
d = 0.005
text(0,d,"$x=0$",ha='center',va='bottom')
text(x[-1],d,"$x=2\\pi$",ha='center',va='bottom')
text(0,-d,"$x_0$",ha='center',va='top')
text(h,-d,"$x_1$",ha='center',va='top')
text(2*h,-d,"$x_2$",ha='center',va='top')
text(x[-3],-d,"$x_{N-2}$",ha='center',va='top')
text(x[-2],-d,"$x_{N-1}$",ha='center',va='top')
xticks([]), yticks([]);

Note that $x_0 = 0$ and the last point $x_{N-1} = 2\pi - h < 2\pi$ ; we do not include the point $x=2\pi$ in the grid since by periodicity, the value there is same as at $x=0$ .

In Python, we can generate the grid points like this.

h = 2 * pi / N
x = h * arange(0,N)

The degree $N/2$ trigonometric polynomial is

I_N u(x) = \sumf \tilu_k \phi_k(x) = \sumf \tilu_k \ee^{\ii k x}

(5)

We will determine the $N$ coefficients $\tilu_k$ from the interpolation conditions

u_j = u(x_j) = I_N u(x_j), \qquad j=0,1,\ldots,N-1

(6)

This is a system of $N \times N$ equations which requires solving a matrix. But we can get the solution without explicitly solving the matrix problem.

First prove the orthogonality relation

\boxed{ \frac{1}{N} \sum_{j=0}^{N-1} \ee^{-\ii p x_j} = \begin{cases} 1, & p = Nm, \quad m=0, \pm 1, \pm 2, \ldots \\ 0, & \textrm{otherwise} \end{cases} }

(7)

Then computing

\begin{align} \sum_{j=0}^{N-1} u(x_j) \ee^{-\ii k x_j} &= \sum_{j=0}^{N-1} \clr{red}{I_N u(x_j)} \ee^{-\ii k x_j} \\ &= \sum_{j=0}^{N-1} \clr{red}{\sum_{l=-N/2}^{N/2-1} \tilde u_l \ee^{\ii l x_j}} \ee^{-\ii k x_j} \\ &= \sum_{l=-N/2}^{N/2-1} \tilde u_l \clr{blue}{\sum_{j=0}^{N-1} \ee^{\ii (l-k)x_j}} \\ &= \sum_{l=-N/2}^{N/2-1} \tilde u_l \clr{blue}{N \delta_{lk}} \end{align}

(8)

we obtain the discrete Fourier coefficients

\boxed{\tilu_k = \frac{1}{N} \sum_{j=0}^{N-1} u_j \ee^{-\ii k x_j}, \qquad k=-N/2, \ldots, N/2 - 1}

(9)

This is known as the discrete Fourier transform (DFT) of $\{ u_j \}$ .

The interpolation conditions

\boxed{u_j = \sumf \tilu_k \ee^{\ii k x_j}, \qquad 0 \le j \le N-1}

(10)

will be called the inverse DFT (IDFT) of $\{ \tilu_k \}$ .

6.12.1.1Lagrange form $^*$ ¶

We can derive a Lagrange form of the interpolation as follows.

\begin{aligned} I_N u(x) &= \sumf \tilu_k \ee^{\ii k x} \\ &= \frac{1}{N} \sumf \sum_{j=0}^{N-1} u(x_j) \ee^{\ii k (x-x_j)} \\ &= \sum_{j=0}^{N-1} u(x_j) \frac{1}{N} \sumf \ee^{\ii k (x-x_j)} \\ &= \sum_{j=0}^{N-1} u(x_j) \psi_j(x) \end{aligned}

(12)

where

\psi_j(x) = \frac{1}{N} \sumf \ee^{\ii k (x-x_j)}

(13)

The functions $\psi_j$ are trigonometric polynomials in $S_N$ that satisfy

\psi_j(x_l) = \delta_{jl}, \qquad 0 \le j,l \le N-1

(14)

Let us define the set of trigonometric functions

S_N = \left\{ v : [0,2\pi] \to \re, \quad v = \sumf a_k \phi_k, \quad a_k \in \complex \right\}

(16)

Define the discrete inner product

\ip{u,v}_N = \frac{2\pi}{N} \sum_{j=0}^{N-1} u(x_j) \conj{v(x_j)}

(18)

and the usual continuous one

\ip{u,v} = \int_0^{2\pi} u(x) \conj{v(x)} \ud x

(19)

In fact, the discrete one can be obtained if we apply Trapezoidal rule to the continuous inner product.

Proof 2

Let

u = \sumf a_k \phi_k, \qquad v = \sumf b_k \phi_k

(21)

Then

\begin{aligned} \ip{u,v} &= \int_0^{2\pi} u(x) \conj{v(x)} \ud x \\ &= \sum_k \sum_l a_k \conj{b_l} \underbrace{\int_0^{2\pi} \phi_k(x) \conj{\phi_l(x)} \ud x}_{2\pi \delta_{kl}} \\ &= 2\pi \sumf a_k \conj{b_k} \end{aligned}

(22)

and

\begin{aligned} \ip{u,v}_N &= \frac{2\pi}{N} \sum_{j=0}^{N-1} u(x_j) \conj{v(x_j)} \\ &= \frac{2\pi}{N} \sum_{j=0}^{N-1} u(x_j) \sumf \conj{b_k} \ee^{-\ii k x_j} \\ &= 2\pi \sumf \left( \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) \ee^{-\ii k x_j} \right) \conj{b_k} \\ &= 2\pi \sumf a_k \conj{b_k} \end{aligned}

(23)

As a consequence, $\ip{\cdot,\cdot}_N$ is an inner product on $S_N$ and

\norm{u}_N := \sqrt{\ip{u,u}_N} = \sqrt{\ip{u,u}} = \norm{u}, \qquad \forall u \in S_N

(24)

Hence the interpolant satisfies

\ip{I_N u - u, v}_N = 0, \qquad \forall v \in S_N

(27)

Thus

$I_N u$ is the orthogonal projection of $u$ onto $S_N$ with respect to the discrete inner product.
The error $u - I_N u$ is perpendicular to $S_N$ and $I_N u$ is the best approximation to $u$ in $S_N$ wrt the discrete norm $\norm{\cdot}_N$ .

This is also related to the property that $\{ \phi_k(x) \}$ is an orthogonal family wrt to the discrete inner product

\ip{\phi_j, \phi_k}_N = 2 \pi \delta_{jk}, \qquad 0 \le j,k \le N

(28)

Proof 4

Starting from the definition of $\tilu_k$

\begin{aligned} \tilu_k &= \frac{1}{N} \sum_{j=0}^{N-1} \clr{red}{u(x_j)} \ee^{-\ii k x_j} \\ &= \frac{1}{N} \sum_{j=0}^{N-1} \clr{red}{ \sum_{l=-\infty}^\infty \hatu_l \ee^{\ii l x_j} } \ee^{-\ii k x_j} \\ &= \sum_{l=-\infty}^\infty \hatu_l \left( \frac{1}{N} \sum_{j=0}^{N-1} \ee^{\ii (l- k)x_j} \right) \end{aligned}

(30)

But by (7)

\frac{1}{N} \sum_{j=0}^{N-1} \ee^{\ii (l-k)x_j} = \begin{cases} 1, & \textrm{if } l-k = Nm, \quad m=0, \pm 1, \pm 2, \ldots \\ 0, & \textrm{otherwise} \end{cases}

(31)

Hence

\tilu_k = \sum_{m=-\infty}^\infty \hatu_{k + Nm} = \hatu_k + \sum_{m=-\infty,m\ne 0}^\infty \hatu_{k + Nm}

(32)

Define

P_N u = \sumf \hatu_k \phi_k \in S_N

(33)

which is a truncation of the Fourier series $Su$ . Also define the summation notation

\sum_{k=-\infty}^\infty = \sumf + \sumfr

(34)

Proof 5

Using Lemma 4, the interpolant can be written as

\begin{aligned} I_N u &= \sumf \tilu_k \phi_k \\ &= \clr{red}{\sumf \hatu_k \phi_k} + \clr{blue}{ \sumf \left( \sum_{m=-\infty,m\ne 0}^\infty \hatu_{k + Nm} \right) \phi_k } \\ &=: \clr{red}{P_N u} + \clr{blue}{R_N u} \end{aligned}

(36)

Hence

u - I_N u = (u - P_N u) - R_N u

(37)

and

u - P_N u = \sumfr \hatu_k \phi_k \qquad \textrm{is orthogonal to} \qquad R_N u \in S_N

(38)

we get

\norm{u - I_N u}_2^2 = \norm{u - P_N u}_2^2 + \norm{R_N u}_2^2

(39)

which proves the desired result.

The result of Lemma 4 and decay properties of Fourier transform are key to showing the error estimate of trigonometric interpolation.

Proof 6

From the smoothness of the function, we have

|\hatu_k| = \order{\frac{1}{|k|^{m+1}}}, \qquad m \ge 1

(41)

The first inequality is easy since

u - P_N u = \sumfr \hatu_k \phi_k \quad\Longrightarrow\quad \norm{u - P_N u}_\infty \le \sumfr |\hatu_k|

(42)

For the second, we use the decomposition

u - I_N u = u - P_N u - R_N u

(43)

so that

\norm{u - I_N u}_\infty \le \norm{u - P_N u}_\infty + \norm{R_N u}_\infty

(44)

But

\norm{R_N u}_\infty \le \sumf \sum_{m=-\infty,m\ne 0}^\infty |\hatu_{k+Nm}| = \sumfr |\hatu_k|

(45)

and hence

\norm{u - I_N u}_\infty \le \sumfr |\hatu_k| + \sumfr |\hatu_k| = 2 \sumfr |\hatu_k|

(46)

The error is given by

\sumfr |\hatu_k| \le C \left[ \frac{1}{(N/2)^{m+1}} + \frac{1}{(N/2+1)^{m+1}} + \ldots \right]

(47)

The sum on the right can be bounded by an integral

\begin{aligned} \frac{1}{(N/2)^{m+1}} + \frac{1}{(N/2+1)^{m+1}} + \ldots &< \int_0^\infty \frac{1}{(N/2 + x - 1)^{m+1}} \ud x \\ &= \frac{1}{m(N/2-1)^{m}} \\ &= \order{\frac{1}{N^{m}}} \quad \textrm{for } N \gg 1 \end{aligned}

(48)

6.12.2Computing the DFT¶

6.12.2.1DFT of a real function¶

If the function $u$ is real, then the DFT comes in complex conjugate pairs

\tilu_{-k} = \conj{\tilu_k}, \qquad 1 \le k \le N/2-1

(49)

and $\tilu_0$ , $\tilu_{-N/2}$ are real. The first property is easy. We have

\tilu_0 = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j)

(50)

which is real, and

\tilu_{-N/2} = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) \ee^{\ii \pi j} = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) \cos(\pi j)

(51)

is also real. When we want to evaluate $I_N u(x)$ at some arbitrary $x$ ,

I_N u(x) = \tilu_0 + \sum_{k=1}^{N/2-1} \left( \tilu_k \phi_k(x) + \conj{\tilu_k \phi_k(x)} \right) + \tilu_{-N/2} \phi_{-N/2}(x)

(52)

the above expression is not real because of the last term. We can modify it as

\begin{align} I_N u(x) &= \tilu_0 + \sum_{k=1}^{N/2-1} \left( \tilu_k \phi_k(x) + \conj{\tilu_k \phi_k(x)} \right) \\ & \qquad + \half \tilu_{-N/2} \phi_{-N/2}(x) + \half \tilu_{-N/2} \phi_{N/2}(x) \end{align}

(53)

i.e.,

I_N u(x) = {\sum_{k=-N/2}^{N/2}}^\prime \hatu_k \phi_k(x), \qquad \hatu_{N/2} = \hatu_{-N/2}

(54)

and the prime denotes that the first and last terms must be multiplied by half. This ensures that the interpolation conditions are still satisfied and $I_N u(x)$ is real for all $x \in [0,2\pi]$ .

However, in actual numerical computation, we may get an imaginary part of order machine precision, and then we have to take the real part. Then, we may as well compute it as before and take real part

I_N u (x) = \textrm{Real} \sumf \tilu_k \phi_k(x)

(55)

6.12.2.2Fast Fourier Transform (FFT)¶

Computing each $\tilu_k$ using (9) requires $\order{N}$ floating point operations (FLOPS). Computing all $N$ coefficients thus requires $\order{N^2}$ FLOPS. This process can be written as a dense $N \times N$ matrix multiplying the vector of values

[u(x_0), \ u(x_1), \ \ldots, \ u(x_{N-1})]^\top

(56)

The matrix has special structure and using some tricks, the cost can be reduced to $\order{N \log N}$ using FFT. It was invented by Gauss and rediscovered 160 years later in Cooley & Tukey, 1965.

6.12.2.3FFT using Numpy¶

Numpy has routines to compute the DFT using FFT, see numpy.fft. It takes the trigonometric polynomial in the form

I_N u(x) = \frac{1}{N} \sum_{k=-N/2+1}^{N/2} \tilu_k \ee^{\ii k x}, \qquad x \in [0,2\pi]

(57)

which is different from our convention.

from numpy import pi,arange
import numpy.fft as fft
h = 2*pi/N
x = h*arange(0,N)
u = f(x)
u_hat = fft(u)

The output u_hat gives the DFT in the following order

\tilu_k = \tilu_0, \ \tilu_1, \ \ldots, \ \tilu_{N/2}, \ \tilu_{-N/2+1}, \ \ldots, \ \tilu_{-2}, \ \tilu_{-1}

(58)

Scipy.fft also provides equivalent routines for FFT.

See more examples of using Fourier interpolation for approximation and differentiation here http://cpraveen.github.io/chebpy.

6.12.3Examples¶

The following function computes the trigonometric interpolant and plots it. It also computes error norm by using ne uniformly spaced points.

# Wave numbers are arranged as k=[0,1,...,N/2,-N/2+1,-N/2,...,-1]
def fourier_interp(N,f,ne=1000,fig=True):
    if mod(N,2) != 0:
        print("N must be even")
        return
    h = 2*pi/N; x = h*arange(0,N);
    v = f(x);
    v_hat = fft(v)
    k = zeros(N)
    n = N//2
    k[0:n+1] = arange(0,n+1)
    k[n+1:] = arange(-n+1,0,1)

    xx = linspace(0.0,2*pi,ne)
    vf = real(exp(1j*outer(xx,k)) @ v_hat) / N
    ve = f(xx)

    # Plot interpolant and exact function
    if fig:
        plot(x,v,'o',xx,vf,xx,ve)
        legend(('Data','Fourier','Exact'))
        xlabel('x'), ylabel('y'), title("N = "+str(N))

    errori = abs(vf-ve).max()
    error2 = sqrt(h*sum((vf-ve)**2))
    print("Error (max,L2) = ",errori,error2)
    return errori, error2

To measure rate of error convergence, we can compute error with $N$ and $2N$ points

E_N = \frac{C}{N^m}, \qquad E_{2N} = \frac{C}{(2N)^m}

(59)

and estimate the rate

m = \frac{\log(E_N / E_{2N})}{\log 2}

(60)

6.12.3.1Infinitely smooth, periodic function $(m = \infty)$ ¶

u(x) = \exp(\sin x), \qquad x \in [0,2\pi]

(61)

f1 = lambda x: exp(sin(x))
fourier_interp(24,f1);

Error (max,L2) =  7.993605777301127e-14 6.463175355583166e-13

6.12.3.2Infinitely smooth, derivative not periodic $(m = 1)$ ¶

u(x) = \sin(x/2), \qquad x \in [0,2\pi]

(62)

g = lambda x: sin(x/2)
e1i,e12 = fourier_interp(24,g)

Error (max,L2) =  0.024834020843920963 0.09773746633443824

e2i,e22 = fourier_interp(48,g,fig=False)
print("Rate  (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))

Error (max,L2) =  0.012408518426188423 0.024996075910098184
Rate  (max,L2) =  1.0009869969530465 1.9672100793199772

6.12.3.3Continuous function $(m = 1)$ ¶

f(x) = \begin{cases} 1 - \frac{2}{\pi} |x - \pi|, & \pi/2 \le x \le 3\pi/2 \\ 0, & \textrm{otherwise} \end{cases}, \qquad x \in [0,2\pi]

(63)

trihat = lambda x: (x >= 0.5*pi)*(x <= 1.5*pi)*(1 - 2*abs(x-pi)/pi)
e1i,e12 = fourier_interp(24,trihat)

Error (max,L2) =  0.03193459816557198 0.15767834270061085

e2i,e22 = fourier_interp(48,trihat,fig=False)
print("Rate  (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))

Error (max,L2) =  0.01586894831194785 0.03976512461360404
Rate  (max,L2) =  1.0089137779952295 1.987408920855071

6.12.3.4Discontinuous function $(m=0)$ ¶

Lemma 6 does not apply in this case.

f(x) = \begin{cases} 1, & \pi/2 < x < 3\pi/2 \\ 0, & \textrm{otherwise} \end{cases}, \qquad x \in [0,2\pi]

(64)

f2 = lambda x: (abs(x-pi) < 0.5*pi)
e1i,e12 = fourier_interp(24,f2)

Error (max,L2) =  0.9957898426758144 2.8589545556697407

e2i,e22 = fourier_interp(48,f2)
print("Rate  (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))

Error (max,L2) =  0.9915466799725409 1.4514563367440128
Rate  (max,L2) =  0.006160606441865679 0.9779865149620691

We do not get convergence in maximum norm but we see convergence in 2-norm at a rate of $1/N$ . There is convergence away from the discontinuities but overall, it is not a good approximation.

References¶

Cooley, J. W., & Tukey, J. W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90), 297–301. 10.1090/S0025-5718-1965-0178586-1
Brunton, S. L., & Kutz, J. N. (2019). Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control (1st ed.). Cambridge University Press. 10.1017/9781108380690

6.12Trigonometric interpolation

6.12.1Interpolation problem¶

6.12.1.1Lagrange form∗^*∗¶