---
title: 11-2. Rectangular spectral collocation
tags: Final reports
---

**Reference**
1. Rectangular spectral collocation, T. A. Driscoll and N. Hale, IMA Journal of Numerical Analysis **36** (2016) 108–132, . [10.1093/imanum/dru062](https://doi.org/10.1093/imanum/dru062)

---

## Questions to be answered

1. Sec3.2: Re-sampling and aliasing

2. Sec4: Rectangular collocation for boundary-value problems

---

### Resampling and aliasing

#### (i) Resampling
Recall the Barycentric interpolation formula :
The polynomial interpolation through $(x_j,f_j)_{j=0}^N$ is given by $p(x) = \frac{\sum\limits_{j=0}^N \frac{\lambda_j f_j}{x-x_j}}{\sum\limits_{j=0}^N \frac{\lambda_j}{x-x_j}}$ with $p(x_j)=f_j,$ where $\lambda_j = (\prod\limits_{k \neq j} x_j-x_k )^{-1}$

In particular, when $\{ x_j \}_{j=0}^N$ are Chebyshev points of second kind, $\lambda_j = \begin{cases}
(-1)^j &,j=1,\cdots ,N-1 \\
\frac{(-1)^j}{2} &,j=0,N
\end{cases}$

From the formula, we can see that the interpolant depends linearly on the values $\{ f_j \}$. That is, the evaluation of $p(y)$ at a set of points $\mathbf{y} = \{ y_i \}_{i=0}^{N-m}$ can be expressed as the matrix-vector multiplication
$$
p(\mathbf{y}) =\begin{bmatrix}
& & \\
& & \\
&P_{N,-m}&\\
& & \\
& & \\
\end{bmatrix} \begin{bmatrix}
f_0\\
f_1\\
\vdots\\
f_{N-1}\\
f_N
\end{bmatrix}
$$
where $P_{N,-m} \in \mathbb{R}^{(N-m+1) \times (N+1)}$ is the **barycentric resampling matrix**
$$[ P_{N,-m} ]_{j,k}=
\begin{cases}
\frac{\lambda_k}{y_j-x_k} (\sum\limits_{l=1}^{N} \frac{\lambda_l}{y_j-x_l})^{-1} &, y_j \neq x_k\\
1 &, y_j = x_k
\end{cases}
$$

**Example :**
$u(x) = e^x \sin 5x$

using LinearAlgebra
using Plots
using FFTW

# resampling f from original second kind chebyshev points to n_new+1 new first/second kind chebyshev points
function Chebresample(n_new, f, kind)
    
    # generate old grid and new grid
    n_old = length(f)-1;
    x_old = -sin.((-n_old:2:n_old)*pi/(2*n_old));
    if kind == 1
        x_new = cos.((1:2:2*n_new+1)*pi/(2*n_new+2));
    else
        x_new = -sin.((-n_new:2:n_new)*pi/(2*n_new));
    end
    
    # construct Barycentric resampling matrix
    P = zeros(n_new+1,n_old+1);
    w = (-1.0).^(0:n_old);  # weight
    w[1] = w[1]/2; w[end] = w[end]/2;
    
    for i = 1:n_new+1
        if x_new[i] in x_old
            P[i,findfirst(x->x==x_new[i],x_old)] = 1;
        else
            temp = sum(w ./(x_new[i]*ones(n_old+1)-x_old));
            for j = 1:n_old+1
                P[i,j] = w[j]/(x_new[i]-x_old[j])/temp;
            end
        end
    end
    
    return x_new, P*f
end

N = 26;
x = -sin.((-N:2:N)*pi/(2*N));
u = exp.(x).*sin.(5*x);

xx=-1:1/100:1;
uu=(exp.(xx).*sin.(5*xx));
NN=13;
(x_new, u_new) = Chebresample(NN,u,1);

plot(xx,uu,label="u(x)")
plot!(x,u,marker=:circle,seriestype = :scatter,label="original sample points",legend=:bottomleft)
plot!(x_new,u_new,marker=:utriangle,seriestype = :scatter,label="resampling points")


u_new_exact = exp.(x_new).*sin.(5*x_new);
print(string("max error = ",norm(u_new_exact-u_new,Inf)));


coefficients :

u_extend = [u; u[N:-1:2]];
u_hat = fft(u_extend)/(2*N);

a = zeros(N+1,1);
a[1] = real(u_hat[1]);
a[2:N] = 2*real(u_hat[2:N]);
a[N+1] = real(u_hat[N+1]);

plot(1:N+1,log10.(abs.(a)),seriestype=:scatter,label="a_k")


#### (ii) Resamlping with differentiation matrix
Recall that $[f_j^{(k)}]=D^k[f_j]$ where $D$ is the differentiation matrix. We can have
$$
p^{(k)}(\mathbf{y})=[P_{N,-m}][f_j^{(k)}]=[P_{N,-m}]D^k[f_j]
$$

**Example :**
$u'(x) = e^x(\sin 5x + 5\cos 5x)$

function cheb(N)
    if N == 0
        D = 0; x = 1;
        return
    end
    x = -sin.((-N:2:N)*pi/(2*N));          # generate grid x_i = cos(i*pi/N)
    c = [2; ones(N-1,1); 2].*(-1).^(0:N);  # c_i*(-1)^i 
    X = repeat(x,1,N+1);                    
    dX = X-X';                             # compute x_i - x_j
    D = (c*(1 ./c)')./(dX+(I(N+1)));       # off-diagonal entries
    D = D - diagm(vec(sum(D,dims=2)));     # diagonal entries
    return D
end

DD = cheb(N);
(x_prime, u_prime) = Chebresample(N-1,D*u,2);

uu_prime = exp.(xx).*(sin.(5*xx)+5*cos.(5*xx));
exact_u_prime = exp.(x_prime).*(sin.(5*x_prime)+5*cos.(5*x_prime));

plot(xx,uu_prime,label="u'(x)")
plot!(x,D*u,marker=:circle,markersize=6,seriestype=:scatter,label="D*u")
plot!(x_prime,u_prime,marker=:utriangle,seriestype=:scatter,color=:yellow,label="resampling points (P*D*u)",legend=:bottomleft)


print(string("max error = ",norm(exact_u_prime-u_prime,Inf)));


#### (iii) Aliasing
Theorem 1 :
For any $N \geq 1$ and $0 \leq k \leq N$, the following Chebyshev polynomials take the same value on the $(N+1)-$point second-kind Chebyshev grid :
$$
T_k,\; T_{2N \pm k},\; T_{4N \pm k},\; T_{6N \pm k},\; \cdots
$$

Theorem 2 :
For any $N \geq 1$ and $0 \leq k \leq N$, the following Chebyshev polynomials take the same value on the $N-$point first-kind Chebyshev grid :
$$
T_k,\; -T_{2N \pm k},\; T_{4N \pm k},\; -T_{6N \pm k},\; \cdots
$$
Now consider downsampling a degree $N$ polynomial $p_{N-1} = \sum\limits_{k=0}^{N-1} a_k T_k(x)$, quoting
a result for second-kind discretizations, to the degree $N-m-1$ polynomial $p_{N-m-1} = \sum\limits_{k=0}^{N-m-1} c_k T_k(x)$ by interpolating at $N−m$ points (where, for simplicity, we restrict $m<\frac{N}{2}$).

It follows from Theorem 1 that
$$
c_k = \begin{cases}
a_k &,\; k=0,1,\cdots,N-2m-2 \text{, and } N-m-1\\
a_k + a_{2(N-m)-(k+1)} &,\; k=N-2m-1,\cdots,N-m-2
\end{cases}
$$ if the $N-m$ points are second-kind Chebyshev points.

And it also follows from Theorem 2 that
$$
c_k = \begin{cases}
a_k &,\; k=0,1,\cdots,N-2m\\
a_k - a_{2(N-m)-(k-1)} &,\; k=N-2m+1,\cdots,N-m-1
\end{cases}
$$ if the $N-m$ points are first-kind Chebyshev points.

Illustration of aliased coefficients : (a) is second to second and (b) is second to first.

**Example :**
$u(x)=e^x\sin5x$

(1) Second-kind to Second-kind

u_extend = [u; u[N:-1:2]];
u_hat = fft(u_extend)/(2*N);

a = zeros(N+1,1);
a[1] = real(u_hat[1]);
a[2:N] = 2*real(u_hat[2:N]);
a[N+1] = real(u_hat[N+1]);

NN=N-1;
(x_new, u_new) = Chebresample(NN,u,2);
u_new_extend = [u_new; u_new[NN:-1:2]];
u_new_hat = fft(u_new_extend)/(2*NN);

c = zeros(NN+1,1);
c[1] = real(u_new_hat[1]);
c[2:NN] = 2*real(u_new_hat[2:NN]);
c[NN+1] = real(u_new_hat[NN+1]);

plot(1:N+1,log10.(abs.(a)),seriestype=:scatter,label="a_k")
plot!(1:NN+1,log10.(abs.(c)),seriestype=:scatter,marker=:utriangle,label="c_k")
plot!(1:NN+1,log10.(abs.(c-a[1:NN+1])),seriestype=:scatter,marker=:+,label="error")


(2) Second-kind to First-kind

u_extend = [u; u[N:-1:2]];
u_hat = fft(u_extend)/(2*N);

a = zeros(N+1,1);
a[1] = real(u_hat[1]);
a[2:N] = 2*real(u_hat[2:N]);
a[N+1] = real(u_hat[N+1]);

NN=N-1;
(x_new, u_new) = Chebresample(NN,u,1);
c = dct(u_new);
c[1] = c[1]/sqrt(NN+1);
c[2:end] = c[2:end]/sqrt((NN+1)/2);

plot(1:N+1,log10.(abs.(a)),seriestype=:scatter,label="a_k")
plot!(1:NN+1,log10.(abs.(c)),seriestype=:scatter,marker=:utriangle,label="c_k")
plot!(1:NN+1,log10.(abs.(c-a[1:NN+1])),seriestype=:scatter,marker=:+,label="error")


From the above illustration, the downsampling of $p_{N−1} =\sum\limits^{N−1}_{
k=0} a_kT_k(x)$ to a grid of $N − 1$ Chebyshev points of the first kind is equivalent to truncating the highest-order Chebyshev coefficient $a_{N−1}$.