1 Interpolation problem ¶ For any even integer N > 0 N>0 N > 0
x j = 2 π j N , j = 0 , 1 , … , N − 1 x_j = \frac{2\pi j}{N}, \qquad j=0,1,\ldots,N-1 x j = N 2 πj , j = 0 , 1 , … , N − 1 The degree N / 2 N/2 N /2
I N u ( x ) = ∑ k = − N / 2 N / 2 − 1 u ~ k ϕ k ( x ) = ∑ k = − N / 2 N / 2 − 1 u ~ k e i k x I_N u(x) = \sumf \tilu_k \phi_k(x) = \sumf \tilu_k \ee^{\ii k x} I N u ( x ) = k = − N /2 ∑ N /2 − 1 u ~ k ϕ k ( x ) = k = − N /2 ∑ N /2 − 1 u ~ k e i k x We will determine the coefficients from the interpolation conditions
u ( x j ) = I N u ( x j ) , j = 0 , 1 , … , N − 1 u(x_j) = I_N u(x_j), \qquad j=0,1,\ldots,N-1 u ( x j ) = I N u ( x j ) , j = 0 , 1 , … , N − 1 Using the orthogonality relation
1 N ∑ j = 0 N − 1 e − i p x j = { 1 , p = N m , m = 0 , ± 1 , ± 2 , … 0 , otherwise \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} N 1 j = 0 ∑ N − 1 e − i p x j = { 1 , 0 , p = N m , m = 0 , ± 1 , ± 2 , … otherwise we obtain the discrete Fourier coefficients
u ~ k = 1 N ∑ j = 0 N − 1 u ( x j ) e − i k x j , k = − N / 2 , … , N / 2 − 1 \tilu_k = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) \ee^{-\ii k x_j}, \qquad k=-N/2, \ldots, N/2 - 1 u ~ k = N 1 j = 0 ∑ N − 1 u ( x j ) e − i k x j , k = − N /2 , … , N /2 − 1 We can derive a Lagrange form of the interpolation as follows.
I N u ( x ) = ∑ k = − N / 2 N / 2 − 1 u ~ k e i k x = 1 N ∑ k = − N / 2 N / 2 − 1 ∑ j = 0 N − 1 u ( x j ) e i k ( x − x j ) = ∑ j = 0 N − 1 u ( x j ) 1 N ∑ k = − N / 2 N / 2 − 1 e i k ( x − x j ) = ∑ j = 0 N − 1 u ( x j ) ψ j ( x ) \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} I N u ( x ) = k = − N /2 ∑ N /2 − 1 u ~ k e i k x = N 1 k = − N /2 ∑ N /2 − 1 j = 0 ∑ N − 1 u ( x j ) e i k ( x − x j ) = j = 0 ∑ N − 1 u ( x j ) N 1 k = − N /2 ∑ N /2 − 1 e i k ( x − x j ) = j = 0 ∑ N − 1 u ( x j ) ψ j ( x ) where
ψ j ( x ) = 1 N ∑ k = − N / 2 N / 2 − 1 e i k ( x − x j ) \psi_j(x) = \frac{1}{N} \sumf \ee^{\ii k (x-x_j)} ψ j ( x ) = N 1 k = − N /2 ∑ N /2 − 1 e i k ( x − x j ) The functions ψ j \psi_j ψ j S N S_N S N
ψ j ( x l ) = δ j l , 0 ≤ j , l ≤ N − 1 \psi_j(x_l) = \delta_{jl}, \qquad 0 \le j,l \le N-1 ψ j ( x l ) = δ j l , 0 ≤ j , l ≤ N − 1 If u ∈ S N u \in S_N u ∈ S N u = ∑ k = − N / 2 N / 2 − 1 c k ϕ k u = \sumf c_k \phi_k u = ∑ k = − N /2 N /2 − 1 c k ϕ k I N u = u I_N u = u I N u = u
Define the discrete inner product
( u , v ) N = 2 π N ∑ j = 0 N − 1 u ( x j ) v ( x j ) ‾ \ip{u,v}_N = \frac{2\pi}{N} \sum_{j=0}^{N-1} u(x_j) \conj{v(x_j)} ( u , v ) N = N 2 π j = 0 ∑ N − 1 u ( x j ) v ( x j ) Let
u = ∑ k = − N / 2 N / 2 − 1 a k ϕ k , v = ∑ k = − N / 2 N / 2 − 1 b k ϕ k u = \sumf a_k \phi_k, \qquad v = \sumf b_k \phi_k u = k = − N /2 ∑ N /2 − 1 a k ϕ k , v = k = − N /2 ∑ N /2 − 1 b k ϕ k Then
( u , v ) = ∫ 0 2 π u ( x ) v ( x ) ‾ d x = ∑ k ∑ l a k b l ‾ ∫ 0 2 π ϕ k ( x ) ϕ l ( x ) ‾ d x ⏟ 2 π δ k l = 2 π ∑ k = − N / 2 N / 2 − 1 a k b j ‾ \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_j}
\end{aligned} ( u , v ) = ∫ 0 2 π u ( x ) v ( x ) d x = k ∑ l ∑ a k b l 2 π δ k l ∫ 0 2 π ϕ k ( x ) ϕ l ( x ) d x = 2 π k = − N /2 ∑ N /2 − 1 a k b j and
( u , v ) N = 2 π N ∑ j = 0 N − 1 u ( x j ) v ( x j ) ‾ = 2 π N ∑ j = 0 N − 1 u ( x j ) ∑ k = − N / 2 N / 2 − 1 b k ‾ e − i k x j = 2 π ∑ k = − N / 2 N / 2 − 1 ( 1 N ∑ j = 0 N − 1 u ( x j ) e − i k x j ) b k ‾ = 2 π ∑ k = − N / 2 N / 2 − 1 a k b k ‾ \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} ( u , v ) N = N 2 π j = 0 ∑ N − 1 u ( x j ) v ( x j ) = N 2 π j = 0 ∑ N − 1 u ( x j ) k = − N /2 ∑ N /2 − 1 b k e − i k x j = 2 π k = − N /2 ∑ N /2 − 1 ( N 1 j = 0 ∑ N − 1 u ( x j ) e − i k x j ) b k = 2 π k = − N /2 ∑ N /2 − 1 a k b k As a consequence, ( ⋅ , ⋅ ) N \ip{\cdot,\cdot}_N ( ⋅ , ⋅ ) N S N S_N S N
∥ u ∥ N : = ( u , u ) N = ( u , u ) = ∥ u ∥ , ∀ u ∈ S N \norm{u}_N := \sqrt{\ip{u,u}_N} = \sqrt{\ip{u,u}} = \norm{u}, \qquad \forall u \in S_N ∥ u ∥ N := ( u , u ) N = ( u , u ) = ∥ u ∥ , ∀ u ∈ S N Hence the interpolant satifies
( I N u − u , v ) N = 0 , ∀ v ∈ S N \ip{I_N u - u, v}_N = 0, \qquad \forall v \in S_N ( I N u − u , v ) N = 0 , ∀ v ∈ S N Thus I N u I_N u I N u u u u S N S_N S N
Using previous Lemma, the interpolant can be written as
I N u = ∑ k = − N / 2 N / 2 − 1 u ~ k ϕ k = ∑ k = − N / 2 N / 2 − 1 u ^ k ϕ k + ∑ k = − N / 2 N / 2 − 1 ( ∑ m = − ∞ , m ≠ 0 ∞ u ^ k + N m ) ϕ k = P N u + R N u \begin{aligned}
I_N u
&= \sumf \tilu_k \phi_k \\
&= \sumf \hatu_k \phi_k + \sumf \left( \sum_{m=-\infty,m\ne 0}^\infty \hatu_{k +
Nm} \right) \phi_k \\
&= P_N u + R_N u
\end{aligned} I N u = k = − N /2 ∑ N /2 − 1 u ~ k ϕ k = k = − N /2 ∑ N /2 − 1 u ^ k ϕ k + k = − N /2 ∑ N /2 − 1 ⎝ ⎛ m = − ∞ , m = 0 ∑ ∞ u ^ k + N m ⎠ ⎞ ϕ k = P N u + R N u Hence
u − I N u = u − P N u − R N u u - I_N u = u - P_N u - R_N u u − I N u = u − P N u − R N u Since
u − P N u = ∑ ∣ k ∣ ≳ N / 2 u ^ k ϕ k u - P_N u = \sumfr \hatu_k \phi_k u − P N u = ∣ k ∣ ≳ N /2 ∑ u ^ k ϕ k is orthogonal to R N u R_N u R N u
∥ u − I N u ∥ 2 = ∥ u − P N u ∥ 2 + ∥ R N u ∥ 2 \norm{u - I_N u}^2 = \norm{u - P_N u}^2 + \norm{R_N u}^2 ∥ u − I N u ∥ 2 = ∥ u − P N u ∥ 2 + ∥ R N u ∥ 2 If u ∈ L 2 ( 0 , 2 π ) u \in L^2(0,2\pi) u ∈ L 2 ( 0 , 2 π ) ∣ u ^ k ∣ = O ( 1 ∣ k ∣ ) |\hatu_k| = \order{\frac{1}{|k|}} ∣ u ^ k ∣ = O ( ∣ k ∣ 1 )
∥ u − P N u ∥ → 0 , ∥ R N u ∥ → 0 \norm{u - P_N u} \to 0, \qquad \norm{R_N u} \to 0 ∥ u − P N u ∥ → 0 , ∥ R N u ∥ → 0 and hence ∥ u − I N u ∥ → 0 \norm{u - I_N u} \to 0 ∥ u − I N u ∥ → 0 N → ∞ N \to \infty N → ∞ u u u
∣ u ( x ) − P N u ( x ) ∣ ≤ ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ → 0 |u(x) - P_N u(x)| \le \sumfr |\hatu_k| \to 0 ∣ u ( x ) − P N u ( x ) ∣ ≤ ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ → 0 For pointwise convergence, we require u ^ k \hatu_k u ^ k
Let m ≥ 2 m \ge 2 m ≥ 2
u u u m − 1 m-1 m − 1 ( 0 , 2 π ) (0,2\pi) ( 0 , 2 π )
u ( m − 1 ) u^{(m-1)} u ( m − 1 )
u ( j ) u^{(j)} u ( j ) 0 ≤ j ≤ m − 2 0 \le j \le m-2 0 ≤ j ≤ m − 2
Then
∥ u − P N u ∥ ∞ ≤ ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ = O ( 1 N m − 1 ) , ∥ u − I N u ∥ ∞ ≤ 2 ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ = O ( 1 N m − 1 ) \norm{u - P_N u}_\infty \le \sumfr |\hatu_k| = \order{\frac{1}{N^{m-1}}}, \qquad
\norm{u - I_N u}_\infty \le 2 \sumfr |\hatu_k| = \order{\frac{1}{N^{m-1}}} ∥ u − P N u ∥ ∞ ≤ ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ = O ( N m − 1 1 ) , ∥ u − I N u ∥ ∞ ≤ 2 ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ = O ( N m − 1 1 ) The first inequality is easy since
u − P N u = ∑ ∣ k ∣ ≳ N / 2 u ^ k ϕ k ⟹ ∥ u − P N u ∥ ∞ ≤ ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ u - P_N u = \sumfr \hatu_k \phi_k \quad\Longrightarrow\quad \norm{u - P_N u}_\infty \le \sumfr |\hatu_k| u − P N u = ∣ k ∣ ≳ N /2 ∑ u ^ k ϕ k ⟹ ∥ u − P N u ∥ ∞ ≤ ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ For the second, we use the decomposition
u − I N u = u − P N u − R N u u - I_N u = u - P_N u - R_N u u − I N u = u − P N u − R N u so that
∥ u − I N u ∥ ∞ ≤ ∥ u − P N u ∥ ∞ + ∥ R N u ∥ ∞ \norm{u - I_N u}_\infty \le \norm{u - P_N u}_\infty + \norm{R_N u}_\infty ∥ u − I N u ∥ ∞ ≤ ∥ u − P N u ∥ ∞ + ∥ R N u ∥ ∞ But
∥ R N u ∥ ∞ ≤ ∑ k = − N / 2 N / 2 − 1 ∑ m = − ∞ , m ≠ 0 ∞ ∣ u ^ k + N m ∣ = ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ \norm{R_N u}_\infty \le \sumf \sum_{m=-\infty,m\ne 0}^\infty |\hatu_{k+Nm}| =
\sumfr |\hatu_k| ∥ R N u ∥ ∞ ≤ k = − N /2 ∑ N /2 − 1 m = − ∞ , m = 0 ∑ ∞ ∣ u ^ k + N m ∣ = ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ and hence
∥ u − I N u ∥ ∞ ≤ ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ + ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ = 2 ∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ \norm{u - I_N u}_\infty \le \sumfr |\hatu_k| + \sumfr |\hatu_k| = 2 \sumfr |\hatu_k| ∥ u − I N u ∥ ∞ ≤ ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ + ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ = 2 ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ The error is given by
∑ ∣ k ∣ ≳ N / 2 ∣ u ^ k ∣ ≤ C [ 1 ( N / 2 ) m + 1 ( N / 2 + 1 ) m + … ] \sumfr |\hatu_k| \le C \left[ \frac{1}{(N/2)^m} + \frac{1}{(N/2+1)^m} + \ldots \right] ∣ k ∣ ≳ N /2 ∑ ∣ u ^ k ∣ ≤ C [ ( N /2 ) m 1 + ( N /2 + 1 ) m 1 + … ] The sum on the right can be bounded
1 ( N / 2 ) m + 1 ( N / 2 + 1 ) m + … < ∫ 0 ∞ 1 ( N / 2 + x − 1 ) m d x = 1 ( m − 1 ) ( N / 2 − 1 ) m − 1 = O ( 1 N m − 1 ) for N ≫ 1 \begin{aligned}
\frac{1}{(N/2)^m} + \frac{1}{(N/2+1)^m} + \ldots &<& \int_0^\infty \frac{1}{(N/2 +
x - 1)^m} \ud x \\
&=& \frac{1}{(m-1)(N/2-1)^{m-1}} \\
&=& \order{\frac{1}{N^{m-1}}} \quad \textrm{for } N \gg 1
\end{aligned} ( N /2 ) m 1 + ( N /2 + 1 ) m 1 + … < = = ∫ 0 ∞ ( N /2 + x − 1 ) m 1 d x ( m − 1 ) ( N /2 − 1 ) m − 1 1 O ( N m − 1 1 ) for N ≫ 1 2 Computing the DFT ¶ If the function u u u
u ^ − k = u ^ k ‾ , 1 ≤ k ≤ N / 2 − 1 \hatu_{-k} = \conj{\hatu_k}, \qquad 1 \le k \le N/2-1 u ^ − k = u ^ k , 1 ≤ k ≤ N /2 − 1 and u ^ 0 \hatu_0 u ^ 0 u ^ − N / 2 \hatu_{-N/2} u ^ − N /2
u ^ 0 = 1 N ∑ j = 0 N − 1 u ( x j ) \hatu_0 = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) u ^ 0 = N 1 j = 0 ∑ N − 1 u ( x j ) which is real, and
u ^ − N / 2 = 1 N ∑ j = 0 N − 1 u ( x j ) e i π j \hatu_{-N/2} = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) \ee^{\ii \pi j} u ^ − N /2 = N 1 j = 0 ∑ N − 1 u ( x j ) e i πj is also real. When we want to evaluate I N u ( x ) I_N u(x) I N u ( x ) x x x
I N u ( x ) = ∑ k = − N / 2 N / 2 ′ u ^ k e i k x I_N u(x) = {\sum_{k=-N/2}^{N/2}}^\prime \hatu_k \ee^{\ii k x} I N u ( x ) = k = − N /2 ∑ N /2 ′ u ^ k e i k x where u ^ N / 2 = u ^ − N / 2 \hatu_{N/2} = \hatu_{-N/2} u ^ N /2 = u ^ − N /2 I N u ( x ) I_N u(x) I N u ( x ) x ∈ [ 0 , 2 π ] x \in [0,2\pi] x ∈ [ 0 , 2 π ]
Numpy has routines to compute the DFT using FFT. It takes the trigonometric polynomial in the form
I N u ( x ) = ∑ k = − N / 2 + 1 N / 2 u ^ k e i k x I_N u(x) = \sum_{k=-N/2+1}^{N/2} \hatu_k \ee^{\ii k x} I N u ( x ) = k = − N /2 + 1 ∑ N /2 u ^ k e i k x 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(v)The output u_hat gives the DFT in the following order
u ^ 0 , u ^ 1 , … , u ^ N / 2 , u ^ − N / 2 + 1 , … , u ^ − 2 , u ^ − 1 \hatu_0, \ \hatu_1, \ \ldots, \ \hatu_{N/2}, \ \hatu_{-N/2+1}, \ \ldots, \ \hatu_{-2}, \ \hatu_{-1} u ^ 0 , u ^ 1 , … , u ^ N /2 , u ^ − N /2 + 1 , … , u ^ − 2 , u ^ − 1 See more examples here http://
3 Examples ¶ The following function computes the trigonometric interpolant and plots it.
# Wave numbers are arranged as k=[0,1,...,N/2,-N/2+1,-N/2,...,-1]
def fourier_interp(N,f,ne=500,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 = int(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(dot(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
f1 = lambda x: exp(sin(x))
fourier_interp(24,f1);Error (max,L2) = 8.01581023779363e-14 4.567409062002914e-13
g = lambda x: sin(x/2)
e1i,e12 = fourier_interp(24,g)Error (max,L2) = 0.024794983262922628 0.0690760697875046
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.012409450043334757 0.017665433009611695
Rate (max,L2) = 0.9986090713583914 1.9672568878077012
def trihat(x):
f = 0*x
for i in range(len(x)):
if x[i] < 0.5*pi or x[i] > 1.5*pi:
f[i] = 0.0
elif x[i] >= 0.5*pi and x[i] <= pi:
f[i] = 2*x[i]/pi - 1
else:
f[i] = 3 - 2*x[i]/pi
return f
e1i,e12 = fourier_interp(24,trihat)Error (max,L2) = 0.0319508570865632 0.11143948223105607
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.015806700636915583 0.02810440673007998
Rate (max,L2) = 1.0153183695992805 1.9873921942542283
f2 = lambda x: (abs(x-pi) < 0.5*pi)
e1i,e12 = fourier_interp(24,f2)Error (max,L2) = 0.9915112319641802 2.0366611449817986
e2i,e22 = fourier_interp(48,f2)
print("Rate (max,L2) = ", log(e1i/e2i)/log(2), log(e12/e22)/log(2))Error (max,L2) = 0.9828401566916884 1.0415601456207677
Rate (max,L2) = 0.0126723113206373 0.9674598167513291
