der2cgl

 DER2CGL    Spectral (Chebyshev Gauss Lobatto) second derivative matrix 

    [d]=der2cgl(x,np) returns the spectral second derivative
    matrix d at the np CGL nodes X (in [-1,1]). np-1 is the polynomial
    degree used. Chebyshev Gauss Lobatto (CGL) grid ORDERED FROM LEFT TO RIGHT


 Input: x = array of CGL nodes on [-1,1] (computed by xwcgl)
        np = number of CGL nodes (=n+1, n=polynomial interpolation degree)

 Output: d = spectral derivative matrix: formula  (2.4.32), pag. 89 CHQZ2

 Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
                    "Spectral Methods. Fundamentals in Single Domains"
                    Springer Verlag, Berlin Heidelberg New York, 2006.

0001 function [d] = der2cgl(x,np) 
0002 % DER2CGL    Spectral (Chebyshev Gauss Lobatto) second derivative matrix
0003 %
0004 %    [d]=der2cgl(x,np) returns the spectral second derivative
0005 %    matrix d at the np CGL nodes X (in [-1,1]). np-1 is the polynomial
0006 %    degree used. Chebyshev Gauss Lobatto (CGL) grid ORDERED FROM LEFT TO RIGHT
0007 %
0008 %
0009 % Input: x = array of CGL nodes on [-1,1] (computed by xwcgl)
0010 %        np = number of CGL nodes (=n+1, n=polynomial interpolation degree)
0011 %
0012 % Output: d = spectral derivative matrix: formula  (2.4.32), pag. 89 CHQZ2
0013 %
0014 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0015 %                    "Spectral Methods. Fundamentals in Single Domains"
0016 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0017 
0018 %   Written by Paola Gervasio
0019 %   $Date: 2007/04/01$
0020 
0021   
0022 
0023 n=np-1; 
0024 x=x(:);
0025 
0026 if n==0 d=0;  return; end; 
0027 
0028 d=zeros(np,np);
0029 dueterzi=2/3;
0030 c=[2;ones(n-1,1);2];
0031 duenp1=2*n*n+1;
0032 d(1,1)=(n^4-1)/15;
0033 for l =2:np
0034 d(1,l)=dueterzi*(-1)^(l-1)/c(l)*(duenp1*(1+x(l))-6)/(1+x(l))^2;
0035 end
0036 for j = 2:n
0037 for l =1:np
0038   if j ~= l          
0039   d(j,l) = (-1)^(j+l-2)/c(l)*(x(j)*(x(j)+x(l))-2)/((1-x(j)^2)*(x(j)-x(l))^2);
0040   else
0041   d(j,j) = -((n*n-1)*(1-x(j)^2)+3)/(3*(1-x(j)^2)^2);
0042   end       
0043 end
0044 end
0045 for l=1:n
0046 d(np,l)=dueterzi*(-1)^(l-1+n)/c(l)*(duenp1*(1-x(l))-6)/(1-x(l))^2;
0047 end
0048 d(np,np)=d(1,1);
0049  
0050 return