How do I form the Chebyshev differentiation matrix in MATLAB?

Question

I have some code that does exactly this, but I do not like to use things I do not understand. Here is the code

N1=N+1;
cheb=cos(pi*(0:N)/N)';
unif=linspace(-1,1,N1)';
if N<3
    x=cheb;
else
    x=cheb+sin(pi*unif)./(4*N);
end
P=zeros(N1,N1);
%N1xN1 zero matrix
xold=2;
%eps= epsilon!
while max(abs(x-xold))>eps

xold=x;   
P(:,1)=1;
P(:,2)=x;
%set first collumn entries to 1 (P(:,1) = 1);
%set second collumn entries to x (P(:,2) = x);
for k=2:N
    P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
    %%Bonnets formula ;)
end
x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1));
end

%---chebyshev differentiation matrix----------------------------------------    ---------
x=-x;
X=repmat(x,1,N1);%set every collumn of X to x
Xdiff=X-X'+eye(N1);
L=repmat(P(:,N1),1,N1);
L(1:(N1+1):N1*N1)=1;
D=(L./(Xdiff.*L'));
D(1:(N1+1):N1*N1)=0;
D(1)=-(N1*N)/4;
D(N1*N1)=(N1*N)/4;

I really want understand, theoretically, how this matrix is formed and then more importantly, implemented in MATLAB.

Up until the comment

% --- chebyshev differenti....

I understand.

I can not seem to find any documentation relating to the formation of this matrix.

NOTE: The syntax I am ok with, it is the actual formation of the matrix I am having trouble with. I have posted this on MSE and it had been put on hold. I had a look on the help centre and this site seemed to fit the description.

Answer 1

I'm assuming that you know how Chebyshev collocation methods work (but if not, let me know and I'll explain a bit more); a good introduction is Nick Trefethen's Spectral Methods in Matlab as well as his Approximation Theory and Approximation Practice (in particular Chapter 21). (But note that this code does Legendre collocation, not Chebyshev collocation!)

Using the Lagrange form of the interpolation polynomial, you can derive an explicit form of the differentiation matrix $D$: $$ D_{ij} = \ell'_j(x_i) = \begin{cases} \frac{\lambda_j}{\lambda_i(x_i-x_j)} & \text{for }i\neq j,\\ \frac{x_j}{1-x_j^2} & \text{for }i= j, \end{cases}$$ where $\ell_j$ are the Langrange nodal polynomials (i.e. $\ell_j(x_k) = \delta_{jk}$) and $\lambda_j$ are the (reciprocal) denominators in the Langrage polynomial $\ell_j$. This is exactly what the final code block computes (in a particularly compact form):

• Xdiff=X-X'+eye(N1) computes all possible differences of $x_i$ and $x_j$ and adds one for $i=j$ (i.e., on the diagonal)
• L=repmat(P(:,N1),1,N1) computes the denominators of $\lambda_j$ from the Chebyshev Vandermonde matrix
• L(1:(N1+1):N1*N1)=1 replaces the diagonal entries with one (since there's no $\lambda_j$ for $i=j$)
• D=(L./(Xdiff.*L')) now computes all entries $D_ij$ simultaneously (it should be obvious for $i\neq j$, the special case $i=j$ is taken care of by the above modifications)
• The remaining lines modify D to account for boundary conditions (the code computes only the matrix corresponding to the interior points).

(For reference, the original code comes from here.)