jacobi_eval

 JACOBI_EVAL: Evaluates Jacobi polynomial P_n^{(\alpha,\beta)} at x\in(-1,1)

              (formula (2.5.3) pag. 92, CHQZ2)
  [p,pd] = jacobi_eval(x,n,alpha,beta) 

 Input: x = scalar or one-dimensional array of length (m) 
        n =  polynomial degree
        alpha, beta= parameters of Jacoby polynomial

 Output: p(m,3) = [P_n^{(\alpha,\beta)}(x),
                   P_(n-1)^{(\alpha,\beta)}(x),
                   P_(n-2)^{(\alpha,\beta)}(x)];

         pd(m,3) = [(P_n^{(\alpha,\beta)})'(x),
                    (P_(n-1)^{(\alpha,\beta)})'(x),
                    (P_(n-2)^{(\alpha,\beta)})'(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 [p,pd] = jacobi_eval(x,n,alpha,beta) 
0002 % JACOBI_EVAL: Evaluates Jacobi polynomial P_n^{(\alpha,\beta)} at x\in(-1,1)
0003 %
0004 %              (formula (2.5.3) pag. 92, CHQZ2)
0005 %  [p,pd] = jacobi_eval(x,n,alpha,beta)
0006 %
0007 % Input: x = scalar or one-dimensional array of length (m)
0008 %        n =  polynomial degree
0009 %        alpha, beta= parameters of Jacoby polynomial
0010 %
0011 % Output: p(m,3) = [P_n^{(\alpha,\beta)}(x),
0012 %                   P_(n-1)^{(\alpha,\beta)}(x),
0013 %                   P_(n-2)^{(\alpha,\beta)}(x)];
0014 %
0015 %         pd(m,3) = [(P_n^{(\alpha,\beta)})'(x),
0016 %                    (P_(n-1)^{(\alpha,\beta)})'(x),
0017 %                    (P_(n-2)^{(\alpha,\beta)})'(x)];
0018 %
0019 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0020 %                    "Spectral Methods. Fundamentals in Single Domains"
0021 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0022 
0023 %   Written by Paola Gervasio
0024 %   $Date: 2007/04/01$
0025 
0026 apb=alpha+beta; ab2=alpha^2-beta^2; 
0027 
0028 nn=length(x);
0029 if nn==1
0030 % x is a scalar
0031 p=1;   pd=0; 
0032 p1=0;  pd1=0;
0033 p2=0;  pd2=0;
0034 if n == 0    
0035    return 
0036 elseif n==1
0037 p1 = p; p2=p1;
0038 p = (alpha-beta+(apb+2)*x)*0.5; 
0039 
0040 pd1 = pd; pd2=pd1;
0041 pd = 0.5*(apb+2); 
0042 else
0043 p1 = p; p2=p1;
0044 p = (alpha-beta+(apb+2)*x)*0.5; 
0045 
0046 pd1 = pd; pd2=pd1;
0047 pd = 0.5*(apb+2); 
0048 for k = 1:n-1 
0049    k1=k+1; k2=k*2; k2ab=k2+alpha+beta;
0050    k2ab1=k2ab+1; k2ab2=k2ab1+1;
0051    p2=p1; p1=p; 
0052    pd2=pd1; pd1=pd; 
0053    a1 = 2*k1*(k1+apb)*k2ab; 
0054 % Gamma(x+n)/Gamma(x) = (x+n-1) * ... * (x+1) * x
0055    a21 = k2ab1*ab2;
0056    a22 = k2ab2*k2ab1*k2ab; 
0057    a3=2*(k+alpha)*(k+beta)*k2ab2;
0058    p = ((a21+a22*x)*p1-a3*p2)/a1; 
0059    pd= (a22*p1+(a21+a22*x)*pd1-a3*pd2)/a1; 
0060 end 
0061 end
0062 
0063 else
0064 % x is an array
0065 
0066 [m1,m2]=size(x);
0067 if m1<m2
0068     x=x';
0069 end
0070 m=max(m1,m2);
0071 p=[ones(m,1),zeros(m,1),zeros(m,1)];   pd=zeros(m,3); 
0072 if n == 0    
0073    return 
0074 elseif n==1
0075 p(:,2) = p(:,1); p(:,3)=p(:,2);
0076 p(:,1) = (alpha-beta+(apb+2)*x)*0.5; 
0077 
0078 pd(:,2) = pd(:,1); pd(:,3)=pd(:,2);
0079 pd(:,1) = 0.5*(apb+2); 
0080 else
0081 p(:,2) = p(:,1); p(:,3)=p(:,2);
0082 p(:,1) = (alpha-beta+(apb+2)*x)*0.5; 
0083 
0084 pd(:,2) = pd(:,1); pd(:,3)=pd(:,2);
0085 pd(:,1) = 0.5*(apb+2);     
0086 for k = 1:n-1 
0087     k2=k*2; k2ab=k2+alpha+beta;
0088    p(:,3)=p(:,2); p(:,2)=p(:,1); 
0089    pd(:,3)=pd(:,2); pd(:,2)=pd(:,1); 
0090    a1 = 2*(k+1)*(k+apb+1)*k2ab; 
0091 % Gamma(x+n)/Gamma(x) = (x+n-1) * ... * (x+1) * x
0092    a21 = (k2ab+1)*ab2;
0093    a22 = (k2ab+2)*(k2ab+1)*k2ab; 
0094    a3=2*(k+alpha)*(k+beta)*(k2ab+2);
0095    p(:,1) = ((a21+a22*x).*p(:,2)-a3*p(:,3))/a1; 
0096    pd(:,1)= (a22*p(:,2)+(a21+a22*x).*pd(:,2)-a3*pd(:,3))/a1; 
0097 end 
0098 end
0099 
0100 end
0101 return 
0102  
0103