How do I form the Chebyshev differentiation matrix in MATLAB?

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.

• Welcome to SciComp.SE! The best person to ask would of course be the code's author (Greg von Winckel, who's now at Sandia). Do you know how spectral methods in general and Chebyshev collocation in particular work mathematically? (If not, I recommend Nick Trefethen's book Spectral Methods in Matlab.) – Christian Clason Jul 17 '16 at 11:17
• I am currently working on a project looking at how we can deal with discontinuities when solving PDE's using the spectral method(I am undergrad). I shall take a close look at that text! I am in shock that you know the code's author just by looking at it! – user21030 Jul 17 '16 at 11:24
• You should also look at chebfun, which implements adaptive (piecewise) Chebyshev approximation and collocation for differential equations -- no need to reinvent the wheel! (I happen to know the author, and recognize his style ;)) – Christian Clason Jul 17 '16 at 11:26
• I have asked a separate question and linked this one in it. – user21030 Jul 26 '16 at 11:18

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

• I am currently reading the relevant material. This answer is more than sufficient. Thanks. – user21030 Jul 17 '16 at 12:21
• Hey Christian. I am having trouble with some code further down, I have done the maths for constructing A, but can not see how this translated to code. – user21030 Jul 26 '16 at 9:49
• Ok brilliant! Sorry I was not sure whether you would answer! – user21030 Jul 26 '16 at 10:01