pnleg2

 PNLEG2    Evaluates the second derivative of Legendre polynomial
           of degree n. 

    [p2,p1,p]=pnleg2(x,n)  evaluates (L_n)''(x), (L_n)'(x), L_n(x)
    at the node(s) x, using the three term relation  (2.3.3), pag. 75 CHQZ2.

 Input: x = scalar or column array
        n = degree of Legendre polynomial

 Output: p2 = scalar or column array (same dimension as x)
             with the evaluation of (L_n)''(x)
         p1 = scalar or column array (same dimension as x)
             with the evaluation of (L_n)'(x)
         p  = scalar or column array (same dimension as x)
             with the evaluation of L_n(x)

 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 [p2,p1,p] = pnleg2 (x, n) 
0002 % PNLEG2    Evaluates the second derivative of Legendre polynomial
0003 %           of degree n.
0004 %
0005 %    [p2,p1,p]=pnleg2(x,n)  evaluates (L_n)''(x), (L_n)'(x), L_n(x)
0006 %    at the node(s) x, using the three term relation  (2.3.3), pag. 75 CHQZ2.
0007 %
0008 % Input: x = scalar or column array
0009 %        n = degree of Legendre polynomial
0010 %
0011 % Output: p2 = scalar or column array (same dimension as x)
0012 %             with the evaluation of (L_n)''(x)
0013 %         p1 = scalar or column array (same dimension as x)
0014 %             with the evaluation of (L_n)'(x)
0015 %         p  = scalar or column array (same dimension as x)
0016 %             with the evaluation of L_n(x)
0017 %
0018 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0019 %                    "Spectral Methods. Fundamentals in Single Domains"
0020 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0021 
0022 %   Written by Paola Gervasio
0023 %   $Date: 2007/04/01$
0024 
0025 nn=size(x);
0026 
0027 if nn==1
0028     % x is a scalar
0029 p=0;p1=p;p2=p;
0030 if n<=1
0031 p=0;p1=p;p2=p;
0032 elseif n==2
0033 p=1;p1=p;p2=p;
0034 else
0035 p1=1;  p2=x;  
0036 p11=0; p21=1; 
0037 p12=0; p22=0;
0038 for k =1:n-1    
0039    p3=((2*k+1)*x*p2-k*p1)/(k+1); 
0040    p31=((2*k+1)*(x*p21+p2)-k*p11)/(k+1); 
0041    p32=((2*k+1)*(x*p22+p21*2)-k*p12)/(k+1); 
0042    p1=p2; p2=p3; 
0043    p11=p21; p21=p31; 
0044    p12=p22; p22=p32; 
0045 end 
0046 p2=p32; p1=p31;p=p3;
0047 end
0048 
0049 % x is an array
0050 else
0051 x=x(:);
0052 p=zeros(nn,1);p1=p;p2=p;
0053 if n<=1
0054 p=zeros(nn,1);p1=p;p2=p;
0055 elseif n==2
0056 p=ones(nn,1);p1=p;p2=p;
0057 else
0058 p1=ones(nn,1);  p2=x;  
0059 p11=zeros(nn,1); p21=ones(nn,1); 
0060 p12=zeros(nn,1); p22=zeros(nn,1);
0061 for k =1:n-1    
0062    p3=((2*k+1)*x.*p2-k*p1)/(k+1); 
0063    p31=((2*k+1)*(x.*p21+p2)-k*p11)/(k+1); 
0064    p32=((2*k+1)*(x.*p22+p21*2)-k*p12)/(k+1); 
0065    p1=p2; p2=p3; 
0066    p11=p21; p21=p31; 
0067    p12=p22; p22=p32; 
0068 end 
0069 p2=p32; p1=p31;p=p3;
0070 end
0071 end
0072 return