gausspd.m

function [A,B] = gausspd(name,x,n,p1,p2,p3)
%GAUSSPD Orthogonal polynomials generated by discrete inner products.
%   [A,B] = GAUSSPD(NAME,X,N,P1,P2,...) returns suitable matrices A and B
%   for usage with CLENSHAW for the polynomials with name NAME with
%   parameters P1, P2, ...
%
%   X is a vector with coordinates and N is the maximum order of the
%   polynomials.
%
%   The polynomials are normalized to unit norm for computational
%   reasons. In addition, the measure has also been normalized to unity.
%
%   No error-checking is performed.
%
%   Name          Parameter(s)    Weight
%   -------------------------------------------------------------
%   'charlier'    a               exp(-a) a^x/x!
%   'krawtchouk'  p, N            bin(N,x) p^x (1-p)^(N-x)
%                                 (x = 0..N)
%   'meixner'     c, b            (1-c)^b c^x (b)_x/x!
%   'chebyshev'   N               1/N (x = 0..N-1)
%   'hahn'        a, b, c         gamma(c-a)gamma(c-b)   ...
%                                 /gamma(c)/gamma(c-a-b) ...
%                                 (a)_x (b)_x/(c)_x/x!
%
%   Here bin(N,x) is the binomial coefficients and (a)_x is Pochhammer's
%   symbol. Use GAUSSQD to check the sanity of the parameters P1, P2, ...
%
%   Examples:
%     x = linspace(0,15);
%     n = 5; a = 5;
%     [A,B] = gausspd('charlier',x,n,a);
%     y = clenshaw(A,B);
%     figure, plot(x,y); title('Charlier polynomials')
% 
%     x = linspace(5,20);
%     c = 0.4; b = 3;
%     [A,B] = gausspd('meixner',x,n,c,b);
%     y = clenshaw(A,B);
%     figure, plot(x,y); title('Meixner polynomials')
%
%   See also CLENSHAW, GAUSSQD.

% S. Engblom 2006-11-06 (Revision)
% S. Engblom 2006-01-12

switch name
 case 'charlier'
  x = reshape(x,1,[]);
  nn = reshape(1:n,n,1);
  sn = sqrt(nn);
  sa = 1./(sqrt(p1)*sn);
  A = [ones(size(x)); ...
       tsum(sa.*(nn+(p1-1)),-sa*x,[1],[1 2])];
  B = -[0; [0; sn(1:end-1)]./sn];
 case 'krawtchouk'
  [A,B] = gausspd('meixner',x,n,-p1/(1-p1),-p2);
 case 'meixner'
  x = reshape(x,1,[]);
  nn = reshape(1:n,n,1);
  nb = sqrt((nn+p2).*(nn+1));
  sc = sqrt(p1);
  A = [ones(size(x)); ...
       tsum((p1*(n+p2)+n)./(sc*nb),(p1-1)/sc./nb*x,[1],[1 2])];
  B = -[0; [0; nb(1:end-1)]./nb];
 case 'chebyshev'
  [A,B] = gausspd('hahn',x,n,1-p1,1,1-p1);
 case 'hahn'
  x = reshape(x,1,[]);
  nn = reshape(0:n,n+1,1);

  w = p3-p1-p2;
  ww = w-2*nn;
  w1 = ww-1;
  w2 = ww+1;

  An = -(ww.*w1)./((w-nn).*(p1+nn).*(p2+nn));
  % *** hopefully, the only singular special case...
  if w == -1
    num = nn.*(w+p1-nn).*(w+p2-nn).*w1;
    den = (w-nn).*(p1+nn).*(p2+nn).*w2;
    Cn = [0; num(2:end)./den(2:end)];
  else
    Cn = (nn.*(w+p1-nn).*(w+p2-nn).*w1)./((w-nn).*(p1+nn).*(p2+nn).*w2);
  end
  A = [ones(size(x)); ...
       tsum(1+Cn(1:end-1),An(1:end-1)*x,[1],[1 2])];
  B = -[0; Cn(1:end-1)];
  % *** implementation of normalization is sub-optimal
  nrm = realsqrt([1; cumprod(Cn(2:end).*An(1:end-1)./An(2:end))]);
  A = tprod(A,[1; nrm(1:end-1)./nrm(2:end)],[1 2],[1]);
  B = B.*[1; 1; nrm(1:end-2)./nrm(3:end)];
 
 otherwise
  error('Unknown polynomial.');
end